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We study the transverse momentum dependent (TMD) evolution of the Collins azimuthal asym¬ 
metries in e + e _ annihilations and semi-inclusive hadron production in deep inelastic scattering 
(SIDIS) processes. All the relevant coefficients are calculated up to the next-to-leading logarith¬ 
mic (NLL) order accuracy. By applying the TMD evolution at the approximate NLL order in the 
Collins-Soper-Sterman (CSS) formalism, we extract transversity distributions for u and d quarks and 
Collins fragmentation functions from current experimental data by a global analysis of the Collins 
asymmetries in back-to-back di-hadron productions in e + e“ annihilations measured by BELLE and 
BABAR Collaborations and SIDIS data from HERMES, COMPASS, and JLab HALL A experi¬ 
ments. The impact of the evolution effects and the relevant theoretical uncertainties are discussed. 

We further discuss the TMD interpretation for our results, and illustrate the unpolarized quark 
distribution, transversity distribution, unpolarized quark fragmentation and Collins fragmentation 
functions depending on the transverse momentum and the hard momentum scale. We make detailed 
predictions for future experiments and discuss their impact. 


I. INTRODUCTION 

The transversity distribution function is one of the three leading-twist quark distributions of nucleon that describe 
completely spin-1/2 nucleon M- Different from the other two, unpolarized and helicity distributions, the quark 
transversity is difficult to measure in experiment because of its chiral-odd nature @ ■ In order to study it in a physical 
process, one has to couple it to another chiral-odd function. The first moments (integral over the longitudinal 
momentum fraction) of the quark transversity distributions lead to the quark contributions to the nucleon tensor 
charge, which is a fundamental property of the nucleon. 

An important channel to investigate the quark transversity distribution is to measure the Collins azimuthal spin 
asymmetries in semi-inclusive hadron production in deep inelastic scattering (SIDIS) Q. Measurements have been 
made by the HERMES Collaboration [(J 0], the COMPASS Colaboration [8j], and JLab HALL A [§| experiments. 
However, the extraction of the quark transversity distributions requires the knowledge of the Collins fragmentation 
functions, which are different from the usual unpolarized fragmentation functions. It was further suggested to measure 
the Collins fragmentation functions from the azimuthal angular asymmetries of two back-to-back hadron productions 
in e + e _ annihilations 0. Recently both BELLE and BABAR Collaborations have studied these asymmetries at 
the B-factories at center of mass energy around yfs ~ 10.6 GeV flll - [l3| . Thanks to the universality of the Collins 
fragmentation functions 0, we will be able to combine the analysis of these two processes to constrain the quark 
transversity distributions. The effort to extract the transversity distributions and Collins fragmentation functions has 
been carried out by the Torino-Cagliari-JLab group extensively in the last few years (l5l4l7| . Transversity coupled to 
the so-called dihadron interference fragmentation functions is employed to study transversity in its collinear version 
in Ref. (l8| . These results have demonstrated the powerful capability of the Collins asymmetry measurements in 
constraining the quark transversity distributions and hence the nucleon tensor charge in high energy scattering exper¬ 
iments. In this study we will implement, for the first time, the appropriate QCD evolution for the phenomenological 
studies of Refs. [ lj| 1 171 and thus improve significantly our understanding of transversity distribution and Collins frag¬ 
mentation functions. We will also show the consistency with previous phenomenological results. A brief summary of 
our results has been published in Ref. 0. 

The appropriate QCD evolution for these low transverse momentum hard processes is the so-called transverse 
momentum dependent (TMD) evolution, which follows from factorization theorems and has been well developed in 
recent years, follo wing the pioneering works by Collins-Soper-Sterman (CSS) 0, HT|. In particular, the Collins 2011 
formalism of Ref. 0] constructs a well defined universal TMDs that “absorb” effects of soft gluon radiation which 
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was traditionally signed out in a separate factor in Refs. [HSU, and defines a hard part function that contains the 
process dependence. This allows for an explicit universality of the TMDs in the phenomenological applications of the 
hard scattering processes mentioned above. 

The TMD evolution effects in the Collins asymmetries in the above processes have been estimated in Refs. [23f . 
The TMD factorization is an important step to derive the results for the physical observables we are interested in, 
and has been shown to be valid for processes with two separate measured momentum scales Qt -C Q, such as SIDIS, 
Drell-Yan and e + e _ annihilation into back-to-back hadrons. Here the small scale Qt corresponds to the measured 
transverse momentum of, for instance, produced hadron in SIDIS, lepton pair in Drell-Yan. TMD factorization is 
formulated in such a way that one can calculate cross sections up to the values Qt ~ Aqcd and thus incorporates 
non-perturbative information on the hadron structure. Schematically the TMD factorization expresses the transverse- 
momentum-dependent differential cross section as a convolution of a so-called hard part H , which is specific for the 
process and thus process dependent, and universal TMD parton distributions and/or TMD parton fragmentation 
functions, collectively called TMDs. These TMDs are universal (for the “naively time reversal odd” functions such as 
Sivers function [24], [25[ and Boer-Mulders function @ the universality is generalized [23, HH) and can be associated 
with properties of specific hadrons. In this sense TMDs represent the three-dimensional partonic structure of the 
incoming nucleons as well as outgoing hadrons. Evolution equations are used to calculate the dependence of TMDs 
on the hard scale Q. Since the definition of TMDs contains the so-called light-cone singularity [53, the detailed 
calculations depend on the scheme to regulate this singularity [2oM22l . l29l , which leads to the scheme dependence 
in the TMD factorization. Although there are different ways to formulate the TMD factorization and to define 
the TMD distribution and fragmentation functions, the energy evolution (historically called “resununation”) for the 
physical observables (including the transverse momentum dependent differential cross sections and spin asymmetries) 
will take the same form in all schemes. In particular, after solving the evolution equations, the final results are 
identical to each other in all TMD factorization schemes, where the TMDs are expressed in terms of their collinear 
counterparts with perturbatively calculable coefficients, and the evolution effects are included in the exponential factor 
- the so-called Sudakov-like form factors. Therefore, in terms of a phenomenological study, all TMD factorization 
and evolution calculations will be identical to that originally proven in the form of CSS [2l| . Interpretation of results 
and individual functions depends of course on the scheme and one should be very careful when giving interpretations. 

TMD evolution is performed in coordinate 6-space, where b is conjugate to the k± in momentum space through the 
Fourier transformation and corresponds to the transverse distance separating the quark/gluon fields. The usage of 
6-space highly simplifies the expressions for the cross sections which become simple products of 6 dependent TMDs 
in contrast to convolutions in k± space. In order to calculate the measured cross-sections (and individual TMDs) one 
performs a two dimensional Fourier transform to the physical Qt (or k±) space. A very unique feature of TMD/CSS 
formalism is the fact that the evolution kernel becomes non-perturbative at large separation distances 6; while at 
small 6 -C 1 /Aqcd it is perturbative and can be calculated order by order in strong coupling constant a s (l/6). 
Over short transverse distance scales, 1/6 becomes a legitimate hard scale, and the 6 dependence of TMDs can be 
calculated in perturbation theory and related to their collinear counterparts, such as collinear parton distribution 
(PDFs), fragmentation functions (FFs) or multiparton correlation functions. The important non-perturabrive part of 
the so-called Soft factor that corresponds to vacuum expectation value of Wilson loops is predicted [53 to be process 
independent, it is also independent of the fact that the individual TMDs is distribution or fragmentation function 
and independent of the particular value of momentum fractions Xb or Zh measured. It may depend on the parton 
type, quark or a gluon, in this paper we are going to consider only quark distribution and fragmentation TMDs. The 
information on intrinsic non-perturbative motion of patrons associated with the hadron wave function is encoded in 
non-perturbative inputs for TMD PDFs and FFs and in turn universal in different processes but in principle dependent 
on the parton/hadron type and on value of Xb or Zh- 

The implementation of the TMD formalism requires parametrization of the non-perturbative inputs [36j - l42| for 
the TMDs. The growth of a s (l/6) at large values of 6 can be tamed by the so-called 6*-prescription (which we 
will follow in this paper) originally introduced in the CSS formalism [21] that allows one to avoid Landau pole in 
strong coupling constant and provides a smooth transition from perturbative to non-perturbative regimes. Fits of 
experimental data utilizing 6* prescription have been well developed in the literature, in particular, in the publications 
of the BLNY -typ e of parameterizations S3- Other choices have been made in the literature, see, for example, 
Refs. [‘171 .1 11 1 71. However, in all these implementations of the TMDs in the CSS formalism, an important step is 
to verify that they provide a robust method of treating non-perturbative physics and can well describe the existing 
experimental data [51, S3 • 

For the Collins asymmetries studied in this paper, we extend the CSS formalism to the azimuthal angular asym¬ 
metries and in the relevant hard processes. This involves the Collins-Soper (CS) evolution equation for the fc^-odd 
distribution and fragmentation functions, which were derived in Refs. [39l. I l8l 1 5.'ll ]. In our calculations, we apply the 
TMD evolution at the approximate next-to-leading-logarithmic order (NLL’) as specified below. The formalism fol¬ 
lows the CSS procedure for the unpolarized cross section, and is similar to that of Sivers asymmetries in SIDIS and 
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Drell-Yan processes [39M4ll. [54| . Wc will derive the perturbative coefficients at one-loop order as well. 

There exists a freedom (scheme-dependence) to separate out the so-called hard factor from the splitting function 
contribution in the CSS formalism [451 ] . This provides a useful way of interpretation for the final results in terms of 
the TMDs [22|, [55], [5f| . It allows to interpret a part of the splitting functions in CSS as a universal TMDs splitting 
functions, and the difference in the coefficients can be regarded as part of hard factors. Once rigorously defined, we 
shall have a unique interpretation of the CSS formalism in terms of TMDs. We will elaborate this interpretation in 
details in our paper. 

In applying the CSS evolution at the NLL order, we relate transversity TMD and Collins FF to the collinear quark 
transversity distribution and the collinear twist-3 fragmentation function and include the DGLAP-type scale evolution 
of the latter two collinear distributions. The evolution of transversity distribution is very well known [57l - l59| while 
the evolution of the twist-3 fragmentation function involves multiparton correlation functions mug , as a common 
feature of higher-twist correlation functions. In the following calculations, we will only keep the homogenous terms in 
the splitting kernel, which is an approximation to the complete evolution equation. To differentiate from the complete 
NLL computation, we denote it as NLL' (an approximate NLL). To achieve this precision we include the most recent 
developments from both theory and phenomenology sides [H, GOMeD . |H, [51;, 53( ;60l-[62|. 

The quark transversity distributions are important ingredients for several other spin related asymmetries. For 
example, they contribute to the azimuthal asymmetries of two-hadron fragmentation processes in SIDIS and e + e - 
annihilations [63j, and single inclusive hadron production at large transverse momentum in pp collisions [64l - [66i ]. 
Future RHIC measurements [67jj are going to explore more phenomena related to transverse spin and ultimately to 
the partonic three-dimensional structure of the nucleon. Our results will provide important cross checks and a step 
further toward a global analysis to all these spin asymmetries associated with the quark transversity distributions. 

The rest of the paper is organized as follows. In Section un we review the theoretical framework for the Collins 
azimuthal asymmetries in SIDIS and e + e _ annihilations and derive the associated TMD evolution results and the 
relevant perturbative coefficients. We also reformulate the resummation formalism in an appropriate way to better 
connect to the recently developed TMD formalism in Section IIII FI In Section IIIII we perform the phenomenological 
studies and focus on the global fit of the quark transversity distribution and Collins fragmentation functions from 
the existing experimental data. We make predictions for future experiments and compare our results with previous 
analyses. Finally, we conclude our paper in Section HV1 


II. THE COLLINS AZIMUTHAL ASYMMETRIES IN SIDIS AND e+e' ANNIHILATION 

In this section, we discuss the asymmetries generated by transversity and Collins fragmentation functions in SIDIS 
and e + e _ annihilation. We apply TMD evolution and represent the differential cross sections, spin-dependent and 
spin-independent ones, in a compact form. 


A. Collins Azimuthal Asymmetries in SIDIS 


In the SIDIS, see Fig. [Q a lepton scatters on the nucleon target, and produces a hadron in the final state, 

e(£) + p(P) —> &{£') + h(Ph ) + X , (1) 


by exchanging a virtual photon = £^ — £' with invariant mass Q 2 = — q 2 . We adopt the usual SIDIS variables 


Q 2 


5 ep — (P + £) , Xb — ^ > 


V = 


P-q = Q 2 
P ■ £ XBS e 


Zh = 


P'Ph 
P-q ’ 


( 2 ) 


with S ep = (£ + P ) 2 the center of mass energy square. The differential SIDIS cross section that includes the Collins 
effect, the sin (ifih + cf> s ) modulation, can be written as mm, 


d 5 a(S±) 

dxBdydzhd 2 Ph± 


CTo(a : B ,y,Q 2 ) 


r? i „• /i i j. \ 2(1 y) j-,sin (</>h+<f> s ) 

F uu + sm(</>fj + (f> s ) ^ F ut 


(3) 


where op = 1+ ^ v ^ , and <j> s and (f>h are the azimuthal angles for the nucleon spin and the transverse momentum 

of the outgoing hadron, respectively. Fuu an( l are gpi n _ av eraged and transverse spin-dependent 

structure functions. The latter is related to convolution of transversity distribution and Collins fragmentation function. 
Ellipsis in Eq.® denote other structure functions that we do not consider in this paper. 
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FIG. 1. Semi Inclusive Deep Inelastic Scattering process (SIDIS) in 7 *P center of mass frame. 


The Collins asymmetry measured experimetnally are related to the structure functions as follows: 


A™(^) = 2<sin(^+^)) 


ao(x B ,y,Q 2 ) 2(1 -y) F^ h+M 
vo{xB,y,Q 2 )i + (i-y) 2 F uv 


(4) 


Note that sometimes experimental results (for instance for COMPASS collaboration) are presented by factoring out 
the so-called depolarization factor Dnn- 


Dnn 


2(1 -y) 

1 + (1 - y) 2 ' 


(5) 


Both structure functions, i.e. Fjjt and p^^h+^F ^ ^ e p enc ] on kinematical variables and on the hard scale Q 2 in the 
reaction. It is important to realize that in order to have reliable calculations of corresponding structure functions one 
needs to take into account appropriate scale dependence which is generated by QCD evolution of TMD distribution 
and fragmentation functions. 

Historically the solution of TMD evolution equations [2l| is presented in the fo-space, where in SIDIS b is Fourier 
conjugate variable to Ph±/zh■ The P/i_i_-dependent structure functions can be formulated in terms of the TMD 
factorization, and they can be (omitting %, z dependencies) written as, 

Fuu{Q',Ph±) = — [ (0 , 2 e l ^ h± '^ Zh F uu (Q;b)+Y U u(Q' 1 Ph±) , 

z h J \ Z7T ) 

^collins(Q;^) = 4 / ■S^ eiPh - L ' S / Zh F£>mns(Q'> b ) + Collins(QI P h±) , 

Z h J ( 27r ) 

where F c “ llins is related to the spin-dependent structure function p^F^h+^F as follows: 

sm(<f>h + </>„) p^+^F = e ^Sl 


9± - 2ef e p x 


F p 

collins ’ 


( 8 ) 


( 6 ) 

(7) 


with the unit vector e x defined in Fig.|T| In Eqs. © and (Q, the first TMD term dominates in Ph± <C Q region, and 
the second so-called Y-factor term dominates in the region of Ph± > Q and assures accuracy of the formula in the 
wide region of Ph±- We will neglect the corresponding Y factors as we will consider only the region of low Phj_/z, 
and thus for spin-averaged and transverse spin-dependent structure functions one has 


Fuu{Q'i Ph±) = 2 


F 


isin (<£/i+0 s ) 


UT 


(Q; Ph±) = — 


d 2 b 

(2^ 

d 2 b 

W? 


e iP h ±-b/z hFuu ( Q . b) 


,iPh _ 


' b/zh Ph±Fcoiu ns (Q] b) 


(9) 

( 10 ) 


one notices that while spin independent structure function is a scalar quantity, the spin dependent structure function 
depends on the transverse direction a = 1, 2; see Eqs. (18IT0D . 
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1. Unpolarized structure function with evolution 

The factorization formula for unpolarized structure function Fuu{Q\Ph. l) is well known and has the following 
interpretation (we choose Ji-Ma-Yuan Jll, Q scheme for the moment) in terms of unpolarized distribution and 
fragmentation functions in the b space [48(: 


Fuu(Q;b) = fi( x B,b-, p,(, v)D h/q {z h ,b\ p,£, p)H(Q/p, p)S(b, p; p) , 

<7 


(ii) 


where ff is the unpolarized TMD distribution, D h / q is the unpolarized TMD fragmentation function, £ 2 = 2(v • 
Pa) 2 / v 2 , C 2 = (2v • Ph) 2 /v 2 , and p 2 = (2v ■ v) 2 /v 2 v 2 represent the light-cone singularity regulation parameters. H 
is the hard factor associated with hard scattering and S is the so-called soft function associated with emission of 
soft gluons. Renormalisation group scale p is arbitrary in full QCD, however in truncated perturbative series it is 
chosen to optimize the convergence in such a way that H does not have large logarithmic contributions, log (Q/p), and 
generically p = C\Q with C\ a parameter of order of 1. We will utilize C 1 = 1 and thus u, = Q in our calculations. 
Depending on different schemes; such as Ji-Ma-Yuan [H,[z3, CSS [I3, HU, or Collins-11 [22}, the dependence on these 
parameters will be different. However, the final results for the structure functions arc independent of the schemes, as 
the actual cross-sections do not depend on the auxiliary parameters. Note that historically H factor is absorbed in 
CSS formulation into the definition of Wilson coefficient functions that relate TMDs to the corresponding collinear 
distributions. The final results for the cross-sections are the same in all schemes. However, a slight difference stems 
from the fact that H functions contain a s (p) with renormalization group scale p while coefficient functions, as will 
be explained below, contain a s (pb) with a dynamical scale /i b . At each order of perturbation series, these differences 
are of a higher order in a s . We will dedicate a separate Section [ill FI where we will discuss TMD interpretation of our 
results and give explicit TMD formulas in TMD Collins-11 [22| formulation for all functions and structure functions 
considered in this paper. 

Let us review the definition and the need of different factors. The TMD quark distributions in SIDIS is defined 
through the following matrix: 


M c 


_ p+ [ -ixj-p+ f d 2 b 

J J (2tt) 2 


Ab-k± 


PS 


,0,&)£t(oo;£)A,(oo;0)#*(0) 


PS 


( 12 ) 


with the gauge link 


C v (oo;£) = exp 


-ig J dX v • A(Xv + £)^ . 


(13) 


This gauge link goes to + 00 , indicating that we adopt the definition for the TMD quark distributions for the SIDIS 
process. The unpolarized quark distribution is projected out from the above matrix as, 

M = \ [fi( x > k-L-hpP * 1 + ■ • •] > ( 14 ) 

fi{x,k±) = ^pTr[7 + 7W] . (15) 

However, the above definition of the quark distribution contains soft gluon contribution, which has to be subtracted 
from the naive definition. In addition, there is light-cone singularity if we take the gauge link along the light-front 
direction v with v 2 = 0. The way to regularize this singularity and subtract soft gluon contribution defines the scheme 
for the TMD factorization. 

In the Ji-Ma-Yuan scheme, the gauge link in the TMD definition is chosen to be slightly off-light-cone, n = 
(l“,0 + ,0j_) —> v = (v _ ,v + ,0_l) with v~ v + . Similarly, for the TMD fragmentation function, v was introduced, 
v = (v - , v -1 ", 0 _l) with i~. Because of the additional directions v and v, there are additional invariants: 

C 2 = (2v ■ P) 2 /v 2 , C 2 = (2v • Ph) 2 /v 2 , and p 2 = ( 2v ■ v) 2 /v 2 v 2 . Accordingly, the soft factor is defined as, 

S^(b) = (0|4 cb ,(6)4 h , Q (6)£„ ab (0)£ 5bc (0)|0). (16) 


With soft factor subtraction, the TMD factorization for the unpolarized structure function can be rewritten as 
Fuu(Q; b) = J2 e l /ijmy^b, b; p, (, p)b[ s ^ Y {z h ,b- p, C, p)H^ iy ( Q/p , p) , 


( 17 ) 
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where the subtracted quark distribution and fragmentation functions are defined as 


/ ijmy >) p, (,p) 

DqJMY( Z h,b] P,C,P) 


fi( z h, b; p, C, f>) 

\/S(b,p\p) 
Dq{z h ,b- 1 p 1 £,p) 
\JS(b, p; p) 


(18) 

(19) 


with the soft factor S subtracted from the original TMDs. After solving the evolution equations and expressing 
the TMDs in terms of the integrated parton distributions, the final expressions for TMDs are obtained by setting 
( 2 = ( 2 = pQ 2 ■ Note that in Eqs. (1181191) we understand the square root in the perturbative sense, i.e. for any 
quantity A = 1 + a± a s + ... one has 1 /VA = 1 — 1/2 aia s — ... . 

On the other hand, as explained in Introduction the new Collins-11 approach |22j is an important improvement of 
the original CSS formalism and includes now operator definition of TMDs, and the soft factor subtraction is taken 
to ensure the absence of light-cone singularities in the TMDs and the self-energy divergencies of the soft factors. 
According to this new scheme, the TMD distribution is defined as 


fi 1CC ( X , b; C,Fi p) = fi(x,b;( F ,p)i 


S n ’ v (b) 


S n ’ n (b)S n ’ v (b ) 


( 20 ) 


where = x 2 (2v ■ Pa) 2 /v 2 = 2 (xP^) 2 e 2yn with y n the rapidity cut-off in Collins-11 scheme. Fragmentation 
functions are defined analogously. The unpolarized structure function takes the form 

Fuu(Q-,b)=Y / 4fi JGG ( x B,b;CF,^D J q cc (z h ,b,CD^)H^ c (Q/p) . ( 21 ) 

9 

One can see from Eqs. (11 71211) the formal expression for structure functions take a very simple “parton model”-like 
form. The underlying TMDs, however, unlike the parton model expression are computed with an appropriate QCD 
evolution procedure. In spite of differences in the schemes to define TMDs in Eqs. (I17I21D the final result for the 
structure functions and the cross-sections is scheme independent and reduces to that of the original CSS. One of 
the advantages of the TMD schemes, such as Collins-11 approach [22| or Ji-Ma-Yuan [lij j7lj, is a possibility to 
define process-independent unpolarized TMDs or account for process dependence in fcj_-odd TMDs directly and study 
individual TMD functions. An important process-independent universal non-perturbative contributions [42j can be 
also studied and global TMD fits that include different processes are possible. 

The corresponding TMDs depend on two scales (or (h>) and p, with their dependence encoded in the TMD 
evolution equations. The rapidity evolution with respect to (/ is given by the Collins-Soper (CS) equation [20]: 


d\nf q (x B ,b- (,p) d\i\D q {z h ,bXD\p) , 

—S75— = -«-"MI'.m), 


( 22 ) 


where K(b , p) is the so-called CS kernel [20]. It can be computed perturbatively for small values of b. The dependence 
on the scale p arises from renormalization group equations for f qi D qi and K: 


dK(b , p) 
din p 

dlnf(x B :b-(,p) 

din p 

d\nD q (z hl b, ( D ;p) 
din p 


- 7 K {oi s (p)) , 

7 F (a s {p),(F/p 2 ) , 

7 d{oi s (p),(d/p 2 ) , 


(23) 

(24) 

(25) 


where functions 7 k, 7 f, and 7 d are anomalous dimensions of K, f q , and D qi respectively. Note that the solution of 
these evolution equations does not depend on the scheme to define TMDs. 

The equations and solutions are discussed at length in Refs. I 20 II 22 I [29l El. ITU. Here we will present and discuss the 
final solution. At low values of b <C 1/A qcd, 1 /b becomes a legitimate hard scale. One introduces [2C| an auxiliary 
scale pb = Co/b , with Cq = 2e~ lE and ~ 0.577 the Euler’s constant. The b dependence of TMDs can be computed 
in the perturbative 1/Q -C b <C 1/A qcd region in terms of the collinear parton distribution and fragmentation 
functions. This region corresponds to the transverse momentum which is large compared to hadronic scale, but still 
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small compared to the hard scale (i.e A qcd < ii « Q). That is TMDs in this region are expressed in terms of 
collinear distributions. This sort of relation will be explained later in the paper. Let us mention that usage of such a 
relation helps to obtain a reliable description of the experimental data. 

The energy evolution of TMDs from the scale gj, to the scale Q is encoded in the exponential factor, exp[—5], with 
the so-called Sudakov-like form factor, the perturbative part of which can be written as 


S V ert{Q,b) 


r Q2 dg 2 

A(a s (g)) In % + B(a s (g)) 

hi A 2 

L m J 


(26) 


where A and B coefficients can be expanded as perturbative series A = (a s /ir) n , B = (a s /n) n . 

In our calculations, we will take A^ l \ A^ and B W for the NLL accuracy. Because this part is spin-independent 
as explained in Introduction, these coefficients are the same as those in unpolarized cross sections [2l| and are given 

by Pair 


A™ = C F , AW 


Cf 


( 67 

7T 2 \ 

10 

~Y 

Ca 

V18 

-e) 

9 


T R n f 


B A) 



(27) 


One can see from Eqs. (191101) that in order to reconstruct the measured cross-section one needs to perform the 
Fourier transform over all values of b. The accuracy of the perturbative solution will deteriorate for large values of 
b. In fact a s (fj,b) will hit the so-called Landau pole which is a good indication of presence of the non-perturbative 
physics. Thus one needs to take into account the non-perturbative behavior of TMDs. The original CSS approach [2lj] 
proposed the so-called ^-prescription that introduces a cut-off value b max and allows for a smooth transition from 
perturbative to non-perturbative region and avoids the Landau pole singularity in a s (gb), 


b=> b* =b/ y/l + b 2 lb 2 max , b max < 1/A QCD , 


(28) 


where b max is a parameter in the prescription. From the above definition, 6* is always in the perturbative region 
where b max is normally chosen to be around 1 GeV^ 1 . With the introduction of 6* in the Sudakov form factor, 
the total Sudakov-like form factor can be written as the sum of perturbatively calculable part and non-perturbative 
contribution 


Ssud(Q', b ) => SpertiQ', K) + <Snp(Q; b) , (29) 

where Sjsip(Q;b) is defined as the difference from the original form factor and the perturbative one. This difference 
should vanish as b — > 0, i.e. in the perturbative region, and thus <Snp(Q; b) has the following generic form 

SNp(Q;b) = g 2 (b)\nQ/Q 0 + g 1 (b) . (30) 

The non-perturbative generic functions gi and g\ have very unique interpretations. In particular <72 includes the 
information on large b behavior of the evolution kernel K . This function does not depend on the particular process, 
it does not depend on the scale and has no dependence on momentum fractions xb, z. This contribution should be 
parametrized phenomenologically and an often-used parametrization is 

92 (b) = g 2 b 2 , (31) 

which proved to be very reliable to describe Drell-Yan data and W ± ,Z boson production in the BLNY-type of 
parameterizations [36j. This gaussian-type parametrization suggests that large b region is strongly suppressed 0] 
and in principle can be unreliable to describe data from lower energies which are more sensitive to moderate-to-high 
values of b. Other parametrizations were proposed in Refs. [381] and [43j. For instance that of Ref. fisj has the form: 

92 (b) = 32 In (J^j , (32) 

and allows to describe simultaneously unpolarized multiplicities from SIDIS measurements by HERMES, low energy 
Drell-Yan as well as Z boson production up to LHC energies. In this paper we will follow the parametrization of 
Ref. [43| for g 2 (b). 

The function g\(b) contains information on intrinsic non-perturbative transverse motion of bound partons, in case 
of distribution TMD it depends on the type of hadron and quark flavor as well as potentially on Xb- In case of 
fragmentation TMD it can depend on Zh and type of the hadron produced, and quark flavor. In other words, 
gi(b) is tied to the particular TMD. Parameters in functions g 2 (b) and gi(b) depend on the cut-off value b max in 
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case ^-prescription is used. The non-perturbative factors could be also defined using different prescriptions; such 
as, for example, matching to perturbative form factors of Ref. [73) or using complex 6 plane integration method of 
Ref. [74j | . In this paper we use the standard CSS ^-prescription method that allows us to compare easily with existing 
phenomenology. 

Therefore with the TMD evolution, TMDs can be expressed as |22], HH, HH, 


ft {sub \x B ,b- Q 2 , Q) = e -S WQ,6.)-s& ^? q [a a {Q))C q ^i 8 f[{x B% p b ), 
D[ sub \z hl b ; Q 2 , Q) = e-^ s »°^^- s »r^V q (a s (Q))C^ q ® D h/j (z h , p b ) , 


(33) 

(34) 


where we explicitly embed the scheme dependence of TMDs from Eqs. (|18I19D in the coefficients T q and T> q . Details 
on these functions will be given in Ref. [|56|. In Ji-Ma-Yuan scheme, 


T q = 1 + C F 

Z7T 


v q = 1 + C F 

Ztt 


1 1 1 2 77 9 

lnp- - In p- — -2 

, 1,2 71-2 o 

lnp- - In p- — -2 


(35) 

(36) 


while in Collins-11 scheme T q = 1 + (7(a 2 ) and V q = 1 + (7(a 2 ). The final result for structure function is p independent 
for Ji-Ma-Yuan scheme, so we set p = 1. In Eqs. (I33I34II . ® represents the convolution in the momentum fraction of 
x or z , 


C q ^®fi(x B ,p b ) = Yl f — c <i^ 

i Jx B X 

fliX’Vb), 

(37) 

A _ r 1 „ 

P'ji—q ® Ph/j ( z h , Pb) = ^ ' / P'j-^—q ( 

.7 Jzh Z 

D h/j( Z ^b)- 

(38) 


The same convolutions will be used for transversity and Collins fragmentation functions with appropriate coefficient 
functions later in the paper. The above coefficient functions are 


Cq<— q ' {X 7 Ph) — 8 q 


+ ^ j- (!-x) 


C qi -g(x,p b ) = —T r x (1 - x) 
7r 


C q 'i~ q (z, p b ) — 5 q ' q 


<5(1 -z) + — f^-(l ~z) + Pq^q(z) In 2 

7T \ Z 


Cg^q(z,p b ) = — f^fz + P g <- q (z) Inz 
TT \ Z 

with the usual splitting functions P q -^ q and P g ^ q given by 

1 + z 2 


Pq^q{z) = Cf 


Pg<-q( z ) ~ P'F 


x + 7^(1 — z ) 
(1 - z)+ 2 

1 + (1 - zf 


(39) 

(40) 

(41) 

(42) 

(43) 

(44) 


The (7-functions are chosen to be universal among different TMD schemes, whereas the functions T q and T> q depend on 
the schemes. In Collins-11 schemes, both factors are equal to 1 up to one-loop order. In the Ji-Ma-Yuan scheme, they 
will depend on p. Again, this p dependence in individual TMDs will be cancelled out by the associated p dependence 
in the hard factor H in Eq. (El when we calculate the structure function Fuu{b , Q). 

Substituting the results of Eqs. (1331341) into the factorization formula Eq. ([TT)) . we can write down the structure 
function Frju hr the 6-space as 

Fuu(Q; b ) = e-Sp-tWW-sS? ? ls (Q, b )Fuu^) , (45) 

with non-perturbative form factor decomposed into the distribution and fragmentation contributions, 

5 ™(Q, 6) = S$p(Q, b) + Sj&(Q, b) , 


(46) 
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which should be determined from the global fit to the SIDIS, e + e , and Drell-Yan data. In the standard CSS 
resummaton which we will follow in this paper, together with the hard factor in the TMD factorization of Eq. m 
the functions F q and D q are absorbed into the C-functions by applying the renormalization group equation for the 
running coupling constant in these two factors (45| . With that, we can write down Fjju(b*) as, 

Fuuih) = (c£^ IS) ® (c]f° IS) <SD h/j (z h ,n b j) , ( 47 ) 

q 


where ^ runs over both quark and anti-quark flavors, J{(xb, Pb) and Db/j{zh, pb ) are the usual unpolarized collinear 
parton distribution function and fragmentation function at the scale pb = Co/&*. We emphasize that the above C- 
coefhcients are the same for all TMD schemes if hard factor H and F q , and V q are absorbed in their definition. 
In particular, in the Ji-Ma-Yuan scheme, the p dependence in U of Eq. (fl7j) . F q in Eq. (l33ll . and V q in Eq. (l34l) 
are cancelled out. In the Collins-11 scheme when the hard factor H is absorbed in the definition of C functions, C 
functions become process dependent and equal to those of the standard CSS scheme. The final expressions for C SIDIS 
and (7 sidis do not depend on p, and they are the same in the Collins-11 scheme, which are also the same as those 
used in the CSS literature [ 75, 76| , 




'ql-q' V"! ru/ — u q'q 

Oft 
IT 


5(1 - x) + 


( C F 


(1 — x) — 2C f 6(1 — x) 


C^ s \x,p b ) = ^T R x( 1-x) 


d^\*,l»)=6* q 


5(1 - z) + y (-^-(1 - z) - 2C f 6(1 — z) + P q ^ q (z) In 2 


C^ is \z,p b ) = y (jyz +P g ^ q (z) In zj . 


(48) 

(49) 

(50) 

(51) 


Of course, there is a freedom to have a separate hard factor in Eq. (H5l) . so that the above C-coefficients will be modified 
accordingly, compare to Eqs. (139I40I41I42I) . This is referred to as scheme dependence [dsj in the CSS resummation. 
For the non-perturbative form factors, we will follow the parameterization of [43], 

SS'p D,8 (<3,6) = S 2 In (£) In + (j, + ft) ^ ■ (52) 

where Q q = 2.4 GeV 2 , for the spin-averaged contribution. In the above parameterization, the parameters g q = gi/2 = 
0.106, q? = 0.84, gh = 0.042 (GeV 2 ) have been determined from the analysis of SIDIS and Drell-Yan processes in 
Ref. jTd . In the fit of Ref. [43], it was found that the non-perturbative form factors do not depend on x. We will use 
the non-perturbative factor of Eq. (1521) in this paper. 


2. Collins structure function with evolution 

Now, we turn to the Collins effects contribution to the single transverse spin asymmetry in SIDIS. We start again 
from the factorized TMD expression in the b space [HI [HJ : 

Collin s (Q;fo) = i^£ L “ / ^ fa) (^,6;p,C, At ) J ff(p,Q/ M ) , (53) 

Q 

where h\ is the TMD quark transversity distribution, Hi h / q is the Collins fragmentation function in the b space and 
is defined (omitting scale dependence) as, 

Hih /q {z h ,b) = j d 2 p_ L e- ip± b p‘tHj- h/q (z h ,p± ) . (54) 

Here H^~ h / q (zh,p±) is the quark Collins function as defined in [5l|, which differs by a factor of (—1 /Zh) from the 
so-called “Trento convention” [78 ], 

Hlh/j( Z hiPx) = ~~H 1 f l ^j[Zh,P±)\ r Tientoj 


(55) 
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with p± the transverse component of the hadron with respect to the fragmenting quark momentum. 

The following model independent relation of Collins fragmentation function H^ h ^ q (zh,pi_) and a twist-3 fragmen- 

A / Q\ 

tation function of quark flavor q to hadron h, can be obtained [51[ : 


\Pi 


#%(**) = I d' P± ^Ht h/j (z h , P± ) . 


One often defines the following so-called first moment of Collins fragmentation function 

IPX I 2 


H^Azh) I Trento = / ^Vl 


i h/j y 


2 z\M% 


Hi h/j (Zh, P±) (Trento • 


(56) 


(57) 


We thus find that 


H h/j(zh) = -^MbH^iZh^Trento- 


X(l), 


(58) 


It is straightforward to show that H^j.{zh) can be written as 


H$ j (z h )=n + £ I d A-e ik+ t-\{Tm« + { 


p+OO 


iDl + J' d(-gF a+ (C ) 


m\p h x) 


x <M|^(0)|0) + ft.c.} 


(59) 


where we have chosen the gauge link in Eq. (|59| going to +oo, and F^ is the gluon field strength tensor and we have 
suppressed the gauge links between different fields and other indices for simplicity. Since the Collins function is the 
same under different gauge links [l4j, [79, 80|, we shall obtain the same result if we replace +oo by — oo in the above 
equation. 

The TMD evolution for the quark transversity and Collins fragmentation functions have been derived in the litera¬ 
ture HH (2(1 |H, [H, [H3| • When expressed in terms of the collinear transversity distribution h\{xB) and the twist-three 

~ /o\ 

fragmentation function H y h ( Zh ), they can be written as 


h\ {sub \x Bl b, p-,Q 2 ,Q) = e i Sp ‘ rt (<2’ 6 *) s np (Q ’ b) Ui q {a s {Q)) 5C q ^ q '® h\ (x B ,Hb), 


II 


i7% ±a M P ;Q 2 ,Q)=l^ h 


—ib° 


= -iSpert(Q, 6 ,)-S^(QX) 


Hc{a s (Q)) SCq’i-q ® H^ q ,{z h ,p b ) , 


(60) 

(61) 


where again, the scheme dependence are in the functions ’Hi q (a s (Q)) and 'H c (a s (Q)). They equal to 1 up to one-loop 
order in Collins-11 scheme. The C-coefficient functions are found to be 


bCq^gt (X, Pb ) = 5 q 'q [(5(1 - x) + 0(a 2 )] 

. , AI I _ _L aS 

'q’<-q V~W ■~o,-~Jq'q 


= <5 q’q [<5(1- z) + ^ (p£<_ q (z)\nz) 


where the function PFAz) has the following form, see Eq. (IA61) . 


(62) 

(63) 


P‘(z) = C F 


2 z 3 

ti -p + - z) 

L (1 ~ z )+ 2 


Substituting the above results into the factorization formula, we obtain the final result for F“ )llins as [21l 

' —ib 


F£ aiBa {Q-,b) = 


2 z h 


g —Spert (Q,b*)—S N p co n ins (Q ,b) (5* ) 


with Ecoiiins (5*) given by 

7f;o|lins(^+ ) 


E 


e 2 q (5Cq^i®h\(xB,Pb)) (<5(7. 


(SIDIS) 

j<-q 


Hff.(z h ,p b 


(64) 

48|, [52 

(65) 

( 66 ) 
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The convolutions are defined in Eqs. (1371381) and the relevant coefficient functions up to the first order in a s (compare 
to Eq. (Ol to determine relation to hard factor H) are given by [Hll. l53l. |62|. [76i]: 


SC { q S ™} S) (Z, Hb)(x, Hb) = Sq'q 


5Cif™ ls) (z,v b ) = 5 q/q 


q <-q 


(5(1 — x) H—- (— 2CfS( 1 — x)) , 

7 r J 

<5(1 - z) + ^ (Pq^q{z) In z — 2C F 5(l - z)) 


(67) 

( 68 ) 


where again, the above C-coefficients contain the contributions from the hard factors in the TMD factorization. The 
hard factor is given in Eq. (IA19I) for Ji-Ma-Yuan scheme and in Eq. (1A2HI) for Collins-11 scheme. 

To achieve the evolution at the NLL order, we have to evaluate both transversity h\{xs , fJ-b) and twist-3 fragmenta¬ 
tion function fib) up to the scale /rj, = co/&*- The evolution for the quark transversity is well-known fsil - t&il ]. 

and we will use the leading order result 


d 

d In /r 2 


hl(xB,n) 



dx 


phi 

q^q 


(i) h\{x B /x,n) , 


(69) 


where the splitting kernel 



C F 


25; 

.(!-*)+ 


+ ^U 



(70) 


Note that since gluon transversity distribution for nucleons does not exist Q, the quark transversity h\ does not mix 
with gluons in its evolution and it evolves as a non-singlet quantity. On the other hand, the evolution equation for 
hfi/j was derived in m® and has a more complicated form. However, if we keep only the homogenous term, we 
can write down the evolution equation as ®], 


d 

9 In g 2 




£ l Zh , 


(71) 


where the splitting kernel Pq<_ q of the homogenous term is given in Eqs. (1A6I64)) . and is the same as that for the 
evolution of the quark transversity function, as pointed out in [60]. We will take this approximation in our numerical 
studies below. In order to differentiate from the complete NLL accuracy we will call it NLL' or approximate NLL. 
For the non-perturbative form factors, we follow the parameterizations of [4dj . 


-S , NPc I oiii„s( ( 9> b ) = 92 In In (J^j + (fl, + ' Jh _ 2 ' Jc ^j b 2 , (72) 

where we assume that the quark transversity follows the same parameterization as unpolarized TMD, but introduce 
an additional parameter to constrain the pj_-dependence in the Collins fragmentation. Therefore, g c will be a free 
parameter in the fit. It is also worthwhile to emphasize that the InQ/Qn-dependent part (i.e. < 72 ln( 6 /b*) in our 
formalism above) is universal for all processes in the initial CSS formalism (2l| as well as in the recent TMD 
formalism of [®. The other contributions in non-perturbative Sudakov form factor are Q-independent and can be 
associated with corresponding TMD distribution and fragmentation functions at an initial scale, see e.g. Ref. [2^.i4l|. 

Finally performing Fourier transforms in Eqs. m and m, we obtain the expressions for both spin-averaged and 
spin-dependent structure functions in the transverse momentum space as 


Fuu(Q'i Ph±) = 


r^sin (<£fc+0s) 
UT 


(Q;Ph±) = 


° dbb 
(2tt) 

1 

2.Zh 


Jo(Ph±b/z h ) e-SpertlQ.ftD-SNF ^(Q^Fuuibt,) , 
f°° dbb 2 


(27r) 


Ji(Ph±b/z h ) e- s P" t ( Q ’ f, *>- s ^=°i»-W’^F comns ( 6 „), 


(73) 

(74) 


with Jo.i the usual Bessel functions. 

Let us comment at this point about the usage of relations to collinear distributions in the structure functions F in 
Eqs. (1471661) . One could in principle solve evolution equations starting at a particular scale Q o instead of introducing 
dynamical scale /i& oc 1/6 and try to extract unknown functions, such as Collins fragmentation function or transversity, 
directly from the data without relying on collinear or twist-3 functions. However such a method has certain difficulties, 
both theoretically and phenomenologically. Theoretical difficulty consists in the fact that if one starts from a fixed 
scale Q 0 then the F function will have potentially large logarithms of the type ln(&Qo) which are obviously not present 
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in fib method due to the choice of fib ~ 1/6. Phenomenologically it might also be difficult to model the unique x,z 
and 6 dependence as contained in the collinear function f(x : fib), which further builds in some dependence on the 
collision energy [37]. Presently there are no successful descriptions of experimental observables simultaneously at 
both low and high energies that use the method with fixed starting scale Qq. The method with the fixed starting 
scale can be applied for processes where the measured scale Q is similar to Qo , namely for processes where the most 
important contribution in cross section comes from b ~ 1 /Q ~ 1/Qo- An example of such a description is a fit of 
Sivers functions in Ref. ^cj 5 IHj]. I n our case, the characteristic scales of SIDIS, Q 2 ~ 2.4 GeV 2 , and e + e _ , Q 2 ~ 110 
GeV 2 , are substantially different. It means that the regions of 6 explored are different and one needs to accurately 
take into account correct 6 dependence of TMDs. That is why in this extraction we will use relations to collinear 
distributions, fragmentation functions, and twist-3 functions. 

By applying the CSS formalism, we utilize the well-established framework of the collinear parton distribution and 
fragmentation functions to parameterize the TMDs at the input scale. For the unpolarized case, this is an obvious 
advantage because of the existing global fits for the integrated PDFs. For the Collins fragmentation function case, it 
is also easier to parameterize TMDs in terms of collinear twist-three function, for which the usual DGLAP evolution 
can be applied. Another important point we want to emphasize is that there are DGLAP-type logarithms in the TMD 
formalism when 6 is small. The CSS formalism is the best way to resum these logarithms, by applying the relevant 
scales (fib) in the associated integrated parton distribution and fragmentation functions. This is an important step to 
help the theory convergence in the perturbative calculations. 


B. Collins Azimuthal Asymmetries in e + e 

In this section we present the formulas for the Collins azimuthal asymmetries in back-to-back di-hadron productions 
in e"*"e _ annihilations, 


e~^ T 6 —} hi -f- h 2 -f- X , (75) 

with center of mass energy S = Q 2 = ( P e + + P e -) 2 , and the final state two hadrons with momenta Phi and Ph 2 , 
respectively. We further identify the longitudinal momentum fractions: Zhi = 2\Phi\/Q. Therefore, Zhi represent the 
momentum fractions in the fragmentation functions which describe the fragmentation processes. Ideally, at leading 
order these two hadrons are produced in a back-to-back configuration. However, the gluon radiation and transverse 
momentum dependence in the fragmentation processes will generate non-zero imbalance between the two hadrons. 

To describe the near-back-to-back imbalance between the two hadrons in e + e _ annihilations, the TMD factorization 
can be used to calculate the differential cross sections. In particular, the Collins fragmentation function will lead 
to a cos 2 (j> azimuthal angular asymmetries between these two hadrons. In the literature, there are two proposed 
experimental methods to investigate the Collins effects in this process: (1) one is to define a thrust axis in e + e - 
annihilation and measure the relative azimuthal angular correlation between the two hadrons in the back-to-back 
two jets, which is referred as A 12 asymmetries; ( 2 ) one is to use one hadron as reference to define the azimuthal 
angle of another hadron (in the back-to-back configuration), which is referred as Aq asymmetries. In the former case, 
we will have to measure two azimuthal angles </>i and <j> 2 , and the Collins effects lead to an azimuthal asymmetry 
proportional to cos(</>i + 1 / 2 ), whereas in the latter case only one azimuthal angle 4 >0 is measured and the Collins 
asymmetry appears as cos(2</o)- In the naive TMD factorization (Born level), both asymmetries can be formulated in 
terms of the Collins fragmentation functions for the hadrons. However, only for the second case, we can immediately 
generalize a QCD factorization in terms of the TMDs. For the first case, a certain modification has to be made to have 
a QCD factorization formula. The reason of this complication is that, in order to describe the case of (1), one has to 
define the jet direction, which is beyond the usual situation of the TMD factorization such as the TMD factorization 
in SIDIS and Drell-Yan lepton pair production. 

In this paper, as a first step, we only consider the second case for the Collins asymmetries in e + e _ annihilations. In 
this measurement, see Fig. [2] the transverse momentum dependence is measured for the hadron (hi) relative to the 
direction of hadron ( 6 . 2 ). The total transverse momentum dependence comes from the TMD fragmentation functions 
for hadron hi and hadron h 2 , plus the soft factor generated from the soft gluon radiation. Again, we focus on the low 
transverse momentum region, where TMD factorization is appropriate and reads j86l . 87] 

d 5 a e+ e~^h 1 h 2 +x _ N c Tra 2 m 
dzhidzh 2 <l 2 Ph±d cos 9 2 Q 2 


(l + cos 2 6 ) Z ^ 2 + sin 2 9 cos(2 <j)o)Z ( 


hi h2 
collins 


(76) 


where 9 is the polar angle between the hadron h 2 and the beam of e + e , 4 >0 is defined as the azimuthal angle of 
hadron hi relative to that of hadron h 2 , i.e. of the plane containing hadrons hi and h 2 relative to the plane containing 
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FIG. 2. e + + e —> hi + h 2 + X process in the frame of method (2). 


hadron h 2 and the lepton pair (see Fig. [2]), and Ph _l is the transverse momentum of hadron hi in this frame. We can 
rewrite the contribution corresponding to -^coiiins Eq. (1761) in the following form: 


sin 2 6>cos(2 (/>o)Z^ 2 ns = sin 2 9 (2e“ef - g'f') 


Z 


h\h2 ot/3 
collins 


(77) 


where the unit vector e x represents the transverse direction of the hadron in the hadron frame and is defined in Fig. [2] 
The tensor structure of this term leads to a cos 2^>o azimuthal asymmetries between the two hadrons. 

The structure functions Z^ h2 and Collins h ave the following form , 


Zuu h3 (Q; P h ±) 
^oUin f(.Q;Ph±) 


1 



_0_ e iP h ±-b/z hl Z hih 2 (Q. b) + Yuu{Q . p h±) j 
_0_^-b/ Zhl ^a 0{ Q. b) + y«f iins( Q. p h±) , 


(78) 

(79) 


where the hrst term depends on the TMD fragmentation functions for the two hadrons and dominates in Ph± / Zhi -C Q 
region, and the second term dominates in the region of Ph±/zhi <1 Q- For cos20o asymmetries, we have an additional 
contribution from gluon radiation [86| associated with spin-averaged fragmentation functions. This contribution does 
not depend on the Collins fragmentation function, and is proportional to P£ ± /Q 2 . It will become important at 
relatively large transverse momentum, and should be included in the above Y terms. However, in the following, 
we only consider the low transverse momentum region Ph _1 <C Q , where this contribution is power suppressed as 
compared to the Collins contributions. In addition, in the experimental measurements, the double ratio of the cos 2 (f> 
asymmetries are reported for di-hadron correlations in e + e _ annihilations, where this contribution is cancelled out. 
Therefore, in the following analysis, we will not include neither the contribution from gluon radiation independent of 
Collins FF nor the Y term. 


1. Experimentally measured Collins azimuthal asymmetries in e + e 


Let us now discuss the definitions of the asymmetries associated with Collins fragmentation functions in the actual 
experimental measurements. Collins function generates cos 2</> 0 modulation in the e + e - cross-section, let us rewrite 
Eq. (1761) as follows: 


d 5 c 


+h\h,2 -\-X 


dzhidzh2d 2 Ph±dcos 9 


2 Q 2 


1 (1 T cos 2 9)Z^ h2 ■ R hltl2 (z h i, z h 2 ,9, P h ±) , 


R 


h 1 hn 


(zhi,z h 2,0,Ph±) = 1 + cos(2</> 0 )- 


sin 9 Z , 


collins 


1 + cos 2 9 Zut! 12 

One could also define analogously the P^j_-integrated modulation 


R hlh2 (z h i,z h2 ,9) = 1 + cos(2</>o) 


sin 2 0 JdP h± P h± Z^ 
1 + cos 2 9 J dP h± P h± Zuu h2 


(80) 

(81) 


(82) 
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In order to eliminate false asymmetries BELLE and BABAR consider the ratios of unlike-sign “U” ( 7 r + 7 r“ + 7 r _ 7 T + ) 
over like-sign “L” ( 7 r + 7 r + + 7 r _ 7 r _ ) or charged “C” ( 7 r + 7 r + + 7 r _ 7 r _ + 7 r + 7 r _ + 7 r _ 7 r + ) pion pairs. In our formalism, 
they can be written as follows: 


R u {z h i,z h2 ,d, P h± ) 
R L (z hl ,z h2 ,9 1 P h± ) 
R u {z hl ,z h2 ,9,P hJ _) 
R G (Zhl,Zh2,0,P h ±) 


1 + cos( 2 (/> 0 ) 
1 + cos( 2 (/>o) 


(sin 2 9) 

(1 + cos 2 9) \ 
(sin 2 9) 


yU yL \ 

^collins ^collins \ 


yL ) 
^uu / 


(1 + cos 2 9) V Z% z 


yU yC n 

■^collins ^collins 


z c J 

uu / 


and likewise for the P^j_-integrated modula: 


R U {zhl, Zh 2 , 9) ^ /« , 1 ( S * n ^ ( / dPh±Ph± ^collins _ / dPh±Ph± -^rollins \ 

R L (z h i,z h2 ,0) - l<PoJ (l + cos 2 0)\ fdP h± P h± ZV u fdP h± P h± Z£ u ) 

R u (z hll z h2 ,9) _ . (sin 2 9) / / dP h± P h± Z G nins _ f dP h± P h± Zg m \ 

R c (z h i,z h2 ,9) + l<PoJ (l + cos 2 8 )\ JdP h± P h± ZV u f dP h ±P h ± Z£ u ) 


where the relevant functions are given by 


yU — yl\ 
^ UU ^ uu 


ry'K 7T - *" 
^UU 


ry L, _ ryT\ 
^ UU ^ UU 


yU 

•^collins 


_ y7r"*~7r . yn 7r"^ 

— ^rnllins I ^ r.t 


''collins 5 


''collins 


+ 3L 

_ ^77r”*"7r"*" 

— ^collins 


z c =z u 

uu uu 


yix 7T 

^collins ? 


■'collins 


— 7 
= Zj, 


collins 


■'collins 


Experimentally measured asymmetries A GL and Aq C are then given by 


Ao L (z h i,z h2 ,9,P h _ l) 
Ao C (z h i,z h2 ,9,P h ±) 


Ao L {z hl ,z h2 ,0) 


A GC (zm, Zh2,9) 


( sin2 6 ) ( Filins •^collins A 

(1 + cos 2 9) \ ZU U Z£ u ) ’ 

(sin 2 9) f Z G uins _ Zg llins \ 

(1 + cos 2 9) { Z" u Zg u J ’ 

(sin 2 9) / f dP h± P h± Z G nins _ f dP h± P h± Z G om \ 

(1 + cos 2 9) V / dP h± P h± Z\[ u f dP h± P h± Z£ u ) ’ 

(sin 2 9) / f dP hJ _P h± 

^collins _ I dPh±Ph± Z collins \ 

(1 + cos 2 9) V JdP h± P h± ZU u J dP hJ _P h± Z£ u J ■ 


(83) 

(84) 

(85) 

( 86 ) 

(87) 

( 88 ) 

(89) 

(90) 

(91) 

(92) 


2. Structure functions 

Corresponding structure functions Z^ h2 and Z^ m 2 ns 
b- space, 


in e + e with QCD evolution 

are defined as Fourier transforms of structure functions in 


zt h2 (Q;P h ±) 




d 2 b 

W? 

d 2 b 

W? 


e ip h^l^zt h2 (Q; b ) , 

„iPh±-b/zhl ( opo p/3 n a P\ 
e \^hX^h± - 9± J 


z ^f(Q-,b) 


(93) 

(94) 


where we only keep the term dominant in the low transverse momentum region. According to the TMD factorization, 
we can write down 


z ui h2 (Q; b ) = E e l D Ch(z h i,bXi^) Di u ;>(z h2 ,b,c2-,^ H e u y 


z coit fm) = , 


(95) 

(96) 


where again Q and p are parameters to regulate the light-cone singularities in the TMD fragmentation functions: 
C 2 = (2 Vi ■ Phi) 2 /v 2 and p 2 = (2i>i • v 2 ) 2 /vfv 2 . One-loop calculations can be performed for the above observables, 
and the relevant hard factors shall follow those in SIDIS calculations. In particular, for the Z uu term, the hard factor 
is the same as that for Drell-Yan lepton pair production, which differs from SIDIS. This happens because in e + e _ 
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annihilation, the virtual photon is time-like q 2 > 0, the same as that in Drell-Yan process, whereas in SIDIS, the 
virtual photon is space-like q 2 < 0. Because of spin-independence of hard interaction in perturbative QCD, the hard 
factor ^coiiins will be the same as Z uu , 


jjC e (JMY) /x-. \ TTI 

H collins .H, 


_ rre+e“(JMV)/ 


(M) 


uu (Q;m) — i + ^ 


c F 


in Ji-Ma-Yuan scheme. For Collins-11 scheme, we will obtain similar results. 


In Q— + In p 2 In E — In p 2 + In 2 p + 27T 2 — 4 


(97) 


H. 


c e oiL (JCC) (Q ;») = Hf u e ~( JCC) (Q; p) = l + 


2n 


3 In Qj — In 2 

p 2 


- 7T 2 -; 


(98) 


Following the previous section, we first derive the evolution results for the TMD unpolarized and Collins fragmentation 
functions and substituting the results into the above factorization formulas of Eqs. (1951961) and obtain, 


Zuu h2 {Q]b) = e -Sp,r t (Q,b,)-s^ ( Q,b ) Z^ h2 (b*) , 


Zcll^(Q;b) = 


—ib a \ (—ib^\ _ 


| 


%Zh2 ) 


— Spert(Q,b») —<S'J,Pcollins(Q>^) '7^1^2 / h 'l 
e ^collins \ U *) 


(99) 

( 100 ) 


Again, the energy evolution effects are explicit in exponential factors, and 


Zt h2 (b*)=Y, e t (A 


\Cq ) ® D hl/i (z h i,p b )) (c£_g ] ®D h2/j {z h2 ,Pb 


=E e * ® i.wo) ® ^SE^ 2 ’^)) 


( 101 ) 

( 102 ) 


where the convolution is defined in Eq. 


and the coefficient functions read: 


Cq %4 E.W-) = <5? 


<5(1 - z) + 


/ Cf 




2 (! - 2) + iVgfc) 111 2 + — (ti 2 - 8 ) 5(1 - z) 


C^f\z,p b ) = ^ + P ff <-,(z) In*) , 


W*) = $ 


9 9 


5(! - 2) + ^ ( PgVgW lnz + E E 8 ) ^ _ ^ 


C F , 2 


(103) 

(104) 

(105) 


with the functions P q <- q: P g <^ q , and P q ^ q given in Eqs. (1431) . (l44l) . and (l64l) . respectively. Again, the above C-coefficient 
functions contain the contributions from the associated hard factors in the TMD factorization, similar to the case of 
SIDIS in the last section (see Appendix A for detailed derivations). On the other hand, the non-perturbative form 
factors are parameterized as 


oe T e 

°NP 


ne^e 

*^NP collins 


Q 


(Q, b) = g 2 In ( f- ) In ( -g- ) + ( -g- + ) b 2 , 

(Q. b) = 32 In ( f- ) In ( E- 


hi h2 
9h 9c 9h 9c 


°hl 


''H2 


(106) 

(107) 


where we have utilized the universality arguments for these parameters. Performing Fourier transforms in Eqs. (1931) 
and (l94l) . we have 


Z uu(Q’, p h±) = — 


dbb 


z hi Jo (2tt) 


MPh±b/z hl ) e -Sp-.W.i.)-s&p- W,‘)^(j t) , 


yhih2 
■^collins 


{Q]Ph±) = 


l 


l 


z\ x 4:ZhlZh2 


dbb 3 

W) 


J 2 (P h ±b/z hl ) e -Vrt, 


(108) 

(109) 


with J 2 the associated Bessel function. 
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III. GLOBAL FIT WITH TMD EVOLUTION 
A. Parametrizations 


As we have seen in previous sections, we have two unknown functions to be extracted from experimental data: 

A / Q\ 

collincar transversity distribution h\ and collinear twist-3 fragmentation function H^' q . The QCD evolution of both 

functions is known, x-dependence of h\ and 2 -dependence of at the initial scale Qo should be parametrized. 

In the global fit, we parameterize the quark transversity distributions as 


h\{x,Q 0 ) = N*x a «{l 


xf 


{a q 


b q ) aq+bq 


a q jbq 
lq Oq 


(fi(x,Q 0 ) + gl(x,Q 0 )) 


( 110 ) 


? are the unpolarized CT10 NLO quark 
S9|. In our parametrization we enforce 


at the initial scale Qo, for up and down quarks q = u,d, respectively, where 
distributions [H] and gf are the NLO DSSV quark helicity distributions 

the so-called Soffer positivity bound ;90] of transversity distribution at the initial scale. This bound is known to be 
valid [59) up to NLO order in perturbative QCD. Soffer bound was criticized in Ref. 911 and was predicted to fail so 
it is very interesting to determine phenomenologically if there are signs of such violation in experimental data. Many 
extractions of transversity, for instance those of Refs. 00, indeed show saturation of Soffer bound for d-quark 
transversity. 

In this study, we assume that all the sea quark transversity distributions are negligible. With more data available 
in the future, we hope we can constrain the sea quark as well, in particular, with the Electron-Ion Collider. We leave 
estimates on possible non-zero sea quarks transversity distributions for future publications. 

Similarly, we parameterize the twist-3 Collins fragmentation functions in terms of the unpolarized fragmentation 
functions, 


Hfi(z,Q 0 ) = N c u z a "{ 1 - z)^D w+/u (z,Q 0 ) , (111) 

Hl%(z, Qo) = N c d z a *( 1 - z)P«D„ +/d (z, Q 0 ) , (112) 

which correspond to the favored and unfavored Collins fragmentation functions, respectively. We also utilize the newest 
NLO extraction of fragmentation functions (92) • The new DSS FF set is capable of describing pion multiplicities 
measured by COMPASS and HEMRES collaborations. In fact it is the only set of fragmentation functions that 
accurately describes COMPASS and HERMES data. The quality of the global fit improved from \ 2 /d.o.f. ~ 2.2 for 
previous DSS NLO [93] to x 2 /d.o./. ~ 1.2 for the new NLO fit [92}. Extractions of LO FFs [93} have yielded a much 
less satisfactory description of the available pion data thus NLO sets ought to be used in extractions of TMDs. NLL 

accuracy allows to utilize this set at NLO. We have verified that results presented here are in complete agreement 

with previously published extraction of Ref. 0. 

The rest can be obtained by applying the isospin relations. We also neglect possible difference of favoured/unfavoured 
fragmentation function of u, d and u , d: 

Hi 3 l /U (z, Qo) = Hl 3 } /u (z, Qq) = Q 0 ) = H unf (z, Q„) , (113) 

H {3 l /d {z, Qo) = H (3 } /d (z, Qo) = Hl 3 } /U (z, Qq) = H fav (z, Qo) . (114) 

Strange quark fragmentation deserves an additional attention. Fragmentation of strange quarks to hadrons is different 
from just “unfavored” fragmentation functions, such as D v +/ u and in order to take this into account we will parametrize 
strange quark “unfavored” Collins fragmenation function as 

^(z, Qo) = H® /atS (z, Qo) = N c d z^( 1 - z)f*‘D„ +/arB (z, Qo) . (115) 

We would like to emphasize that in the fit, we will solve the DGLAP evolution equations for both transversity 
and Collins FF to the scale gb = co/&*, in order to be complete at the NLL' order. Numerical solution of DGLAP 
equations is performed in x-space by HOPPET evolution package [94| . Original code of HOPPET is modified by us so that 
transversity splitting functions are included, the initial scale for the evolution is chosen to be Q\ = 2.4 GeV 2 . and the 
HOPPET code is executed using a s (Qo) = 0.327. In our numerical calculations we consistently use 2-loop order result 
for a s {n) with n/ = 5 effective quark flavors and A qcd = 0.225 GeV such that a.AMz) = 0.118. 

For the non-perturbative form factors, we use the following parameters from Ref. j43|: g q = gi /2 = 0.106, 52 = 0.84, 
gh = 0.042 (GeV 2 ). The NLL formula has a large negative contribution coming from C ^ oc — 2Cf8{1— x), see Eq. (1481) . 
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and H W oc —8 Cf, see Eq. (11391) . this makes the need for potentially large K factors in description of the data. We 
assume that K factors will be largest in lowest Q 2 region where a s is relatively large. In fact the fit of SIDIS and 
Drell-Yan data of Ref. [13 revealed Ksidis ~ 2 for COMPASS and HERMES and Kby ~ 1 for Drell-Yan data. In 
asymmetry K factors cancel, so we will not use them in present analysis. 

The existing experimental data does not allow to determine precisely shapes of all polarised distributions in coor¬ 
dinate space, we make a simplifying assumption and allow for the Collins fragmentation function to modify its shape 
with respect to unpolarised fragmentation distributions and have g c as a free parameter in the fit. 

Therefore, we have total of 13 parameters in our global fit: N^, N%, a u , ad, b u , bd, N£, N2, a u , ad, Pd, Pu, 9c 
(GeV 2 ). 


B. Experimental data 

Let us also discuss available experimental data. In this paper we extract Collins fragmentation functions for pions 
and transversity distributions for u, d and favoured/unfavoured Collins fragmentation functions for pions. Thus we 
will select the data involving pion production only. 

The HERMES Collaboration measured Collins asymmetries in electron proton scattering at the laboratory electron 
beam energy 27.5 GeV in production of tt + , tt~, and tt° [95]. The data are presented in bins of xb, Zh, and Ph± 
respectively. Clear non-zero asymmetries were found for both 7r + and tt~ . Large negative asymmetry for tt~ suggest 
that unfavoured Collins fragmentation function is big and not suppressed with respect to the favoured one. 

The COMPASS Collaboration uses muon beam of energy 160 GeV and have measured Collins asymmetries on both 
NH 3 (proton) [96| and LiD (deuterium) [97j targets. The data are presented as function of xb, Zh, and Ph±- Results 
on the proton target are compatible with HERMES findings and asymmetries are found to be compatible with zero 
on the deuteron target. The beam energy of COMPASS is higher than the energy of HERMES and thus COMPASS 
reaches lower values of Xb ~ 10 -3 . For each point in xb the scale Q 2 is higher at COMPASS as one has Q 2 = sxy. 
Both experiments consider Q 2 > 1 GeV 2 in order to be in the perturbative region and the energy of y*p, W 2 > 10 
GeV 2 for HERMES and W 2 > 25 GeV 2 for COMPASS in order to be outside of the resonance region. The COMPASS 
Collaboration considers Zh > 0.2 region and the HERMES Collaboration uses 0.2 < Zh < 0.7 in order to minimize 
both target fragmentation effects and exclusive reaction contribution. All other experimental cuts are described in 
Refs. [MllS- 

Jefferson Lab’s HALL A published data of ^ pion production in 5.9 GeV electron scattering on 3 He (effective 
neutron) target [9:]. Jefferson Lab operates at relatively low energy and reaches higher values of Xb ~ 0.35. 

Information on Collins fragmentation functions is contained in e + e _ at the energy yfs — 10.6 GeV data of the 
BELLE [3 and the BABAR [98j Collaborations. Note that usual feature of TMD evolution is widening of distributions 
with increase of the hard scale. Thus it is very important to check our knowledge against available data on Phx 
distributions and corresponting Ph± dependencies of asymmetries. For this reason, we include BABAR [9S<] data on 
Ph± dependence in our fit. We will also present predictions of Ph± dependence of the unpolarised cross section that 
will be the ultimate test of the model. As we mentioned in Sec. IIIBl we will use Aq data on Collins asymmetries in 
e + e _ in our fit. Both BELLE and BABAR Collaborations require the momentum of the virtual photon Ph±/zhi < 3.5 
GeV in order to remove contributions from hadrons assigned to the wrong hemisphere and it also helps to remove 
contribution from gluon radiation. The analysis of BELLE is performed in (zhi,Zhi) bins with boundaries at Zhi = 
0.2, 0.3, 0.5, 0.7 and 1.0. The BABAR Collaboration chooses 6 z^-bins: [0.15 — 0.2], [0.2 — 0.3], [0.3 —0.4], [0.4 —0.5], 
[0.5 — 0.7], [0.7 — 0.9]. A characteristic feature of the asymmetry is growth with Zhi which is compatible with Zhi 
dependence of theoretical formula, and the asymmetry should vanish in the limit Zhi —> 0. 


C. Fitting procedure 


We proceed with a global fit of SIDIS and e + e data using MINUIT package [991 . 1100 ] by minimizing the total \ 2 : 

N Nt 


X 2 ({a}) = ZT , ma]) ~ E ’ )2 


*=1 1=1 




(116) 


for i = 1,..., N data sets each containing N t data points. Experimental measurement of each point is Ej, experimental 
uncertainty is A Ej and theoretical estimate, Tj, is calculated for a given set of parameters {a}={lV£, N%, a u , ad, b u , 
bd, N£, N'j, a u , ad, Pd, Pu, 9c}- We include both statistical and systematical experimental uncertainty in quadrature. 
Normalization uncertainties are not included in this fit. We have in total N = 26 sets of which Nsidis = 20 sets for 
SIDIS and N e + e - = 6 sets in e + e _ . The formalism is valid for low values of Phj_/z ■C Q, so we include only SIDIS data 
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for P^x-dependence using a conservative choice Ph± < 0.8 GeV. We also limit Ph±/zhi < 3.5 GeV from BELLE fl2| 
and BABAR [98j data following the experimental cuts. The number of points is N^J.^ 13 = 140 and = 122. The 

number of fitted parameters, 13, is adequate for fitting the total number of data points, Ntotai = 262. More flexible 
parametrizations will be explored in future publications. In the fit we use the average values of ( XB),{zh),{y),{Ph±) 
for each bin in SIDIS and {zhi) ,(zh 2 ) ,(Ph±), (sin 2 9) / ( 1 + cos 2 9) for each bin in e + e _ . 

We present an estimate at 90% confidence level (C.L.) interval for t he n ucleon tensor charge contributions and 
estimate errors on our results using the strategy outlined in Refs. jlOll Il02| . The method consists of exploring the 
parameter space {aj} by exploring possible values of y 2 so that 

X 2 (K» <xi.in+ A X 2 , (117) 


where Ay 2 corresponds to the so-called fit tolerance T = \J Ay 2 . In the ideal case of uncorrelated measurements 
without unknown sources of error and gaussian errors of the measured observables, the 68 % C.L. corresponds to 
Ay 2 = 1 and 90% C.L to Ay 2 = 2.71. In typical measurement of asymmetries or other observables one encounters 
either correlated measurements or some inconsistent data sets due to uncontrolled experimental and/or theoretical 
errors. In order to deal with those issues the tolerance is changed with respect to the standard values. 

A very rough idea of a good fit of the data set that contains N points is the resulting y 2 being in the range of 
N±V2N. A more precise quantification of the allowed tolerance or Ay 2 can be estimated by assuming that calculated 
y 2 follows the y 2 -distribution for N degrees of freedom with the probability density function 


1 

2r(A/2) 


£x 2N 

N/ 2-1 

x 2 ’ 

u. 

) exp 

2 


(118) 


The most probable value is the 50% percentile £50 (compare to the goodness of fit): 


r € 50 


d X ' 


2T(A/2) 


/x 2N 

N/2—1 

x 2 ’ 

(t. 

| exp 

2 


= 0.5. 


This percentile is of order of N. The 90% percentile, £ 90 , is accordingly 


rt. 90 


dr 


1 


2T(A/2) 


( x 2 ' 

N/2-1 

x 2 ’ 

(x 

j exp 

2 


= 0.9. 


(119) 


( 120 ) 


The Ay 2 is defined then as 


Ay 2 = £90 — £50 • ( 121 ) 

Analogously we can define 

AX68 = 638 - £50 , ( 122 ) 

for 68 % C.L. In our particular case with 13 fitting parameters we have N = N to tai — 13 = 249, £50 = 248.3, and 
£90 = 278.0, and thus Ay 2 = 29.7. It is comparable to v 2 2 N = 22.3. We also have Ayg g = 10.6. 

For each set of experimental data i, the 90% C.L. is defined as in Ref. 1021 ] 

X? < f £ 90 , ( 123 ) 

note that the value of £90 is renormalized by X' 2 m m/£50 due to the fact that in the total global minimum y 2 min = 
Yli Ximin the va lue of Ximin m ay be away from possible minimal value. 

In order to estimate errors on parameters and on calculation of asymmetries we will utilize a Monte Carlo sampling 
method explained in Ref. 0 - We are going to generate samples of parameters {a^} in the vicinity of the minimum 
found by MINUIT {ao} that defines the minimal value of total Xmin ■ l 11 order to account for correlations in parameters 
and improve numerical performance we will generate correlated parameter samples using CERNLIB’s0 Monte Carlo 


1 http://cernlib.web.cern.ch/cernlib/ 
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generators CORSET and CORGEN utilizing correlation matrix found by MINUIT. We generate 135 sets of parameters 
{cti} that satisfy 

X 2 (K}) <xLin + A X 2 , (124) 

with Ay 2 from Eq. (112111 . By definition these sets correspond to the hyper volume in parameter space that defines 
90% C.L. region. Any observable then will be calculated using these sets and the maximum and minimal value found 
will define our 90% C.L. estimate. This Monte Carlo method underestimates the errors due to the limited number of 
generated parameter sets (135). The errors on asymmetries and functions that we quote are thus estimates and we 
will use a more robust method to estimate the errors on tensor charge. Errors on the tensor charge will be calculated 
using the evaluation of the x 2 -profile by varying parameters of the model and careful analysis of the possible values. 


D. Results 

The resulting parameters after minimization procedure are presented in Table [0 Only relative sign of transversity 
can be determined and we present here a solution with positive u quark transversity as in Refs. jl5l - |T7l l63i . Il03j . 
Indeed transversity and helicity distributions can be related via boost and rotation of correspondding operators, 
however boost and rotation do not commute in quantum theory, thus these two distributions are independent and in 
principle different. It is unlikely that they differ by sign, thus we choose the same sign for u-quark transversity and 
it-quark helicity distribution [89( which is positive. Transversity of d quark is negative. Favored and unfavored Collins 
FFs are of opposite signs, indeed N° < 0, N% > 0 and of approximately the same magnitudes. It means that favored 
Collins fragmentation fu nction is p ositive and unfavored Collins fragmentation function is negative, see Eq. (1581) . The 
corresponding sum rule |l04 ll05| for Collins fragmentation functions read: 

/•I 

EE/ dZf l ZhH 1 fo/j{Zh) |Trento — 0) (125) 

h S h 

which suggests the compensation of favored and unfavored Collins fragmentation functions. 

We observe that parameters that define z dependence of Collins FF a u and a d are different, thus the z-shapes of 
favorite and unfavorite Collins FFs are different. The same is true for transversity distributions, both large-a; region 
controlled by b u and bd and low-a; region controlled by a u and ad indicates that the ir-shape of transversity for u and 
d quarks is different. It might be well possible that fcj_-shape of transversity and Collins fragmentation functions is 
also flavor dependent, however the current experimental data does not allow to determine whether it is true or not. 



= 0.85 ± 0.09 

a u = 

0.69 

± 

0.04 

bu 

= 0.05 

± 

0.04 

fy? = 

= -1.0 ±0.13 

ad = 

1.79 

± 

0.32 

b d 

= 7.00 

± 

2.65 


= -0.262 ±0.025 

CX-u = 

1.69 

± 

0.01 

fdu = 

= 0.00 

± 

0.54 

N c d -- 

= 0.195 ±0.007 

O'd = 

0.32 

± 

0.04 

Pd 

= 0.00 

± 

0.79 

g c = 

= 0.0236 ± 0.0007 

(GeV : 

! ) 








Xmin = 218.407 Xmin/n.d.O.f = 0.88 

TABLE I. Fitted parameters of the transversity quark distributions for u and d and Collins fragmentation functions. The fit 
is performed by using MINUIT minimization package. Quoted errors correspond to MINUIT estimate. 


The total = 218.407 and x 2 / n d.o.f. = 218.407/249 = 0.88. We calculate the goodness of fit using the well 
known formula [1061 ]: 


P(Xmim n d.o.f.) — 1 


d X 2 


2r(n(i. 0 ././2) 


nd.o.f. / 2—1 


exp 


21 


(126) 


The goodness of fit describes how well it fits a set of observables. In principle if the model adequately describes the 
data then one would expect x 2 / n d.o.f. — 1- In case x 2 / n d.o.f. 1 the model fails to describe the data, x 2 / n d.o.f. "C 1 
means that the model starts fitting the statistical noise in the data. Notice that in our fit we obtained x 2 / n d.o.f. = 0.88 
which means that the number of parameters is adequately chosen. An attempt to extract more information from the 
data, such as flavor dependence etc, would lead to x / n d,.o.f. ■C 1. One of course can estimate a number of different 
hypothesis and we leave those estimates for further publications. 
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We obtain that the probability that the fit indeed is the underlying mechanism for the measured asymmetries is 
P(218.407,249) = 92%. This gives us full confidence in presented results. It is very important to note that we have 
used the data from two different processes implementing appropriate factorization and evolution. Thus we have also 
presented a phenomenological proof that these processes, SIDIS and e + e _ , are consistent with TMD factorization. 

The results of the fit including partial values of y 2 are presented in Table [TT] for SIDIS experiments and in Ta¬ 
ble IIIII for e + e _ experiments. One observes that y 2 values are quite satisfactory and homogeneous for both SIDIS, 
Xsidis/ ^toiM S = 0-93, and e + e _ , x1+ e - /^totai = 0.72. TMD factorization approach is describing data of both 
SIDIS and e + e _ adequately. 

Description of SIDIS data is very good. The data span the energy range starting from P^ AB = 5.9 to P;™ RMES = 
27.5, and to -F)oP MPASS = 160 GeV. The resolution scale changes also in a relatively wide region 1 < ( Q 2 ) < 6 (GeV 2 ) 
for HERMES and 1 < ( Q 2 ) < 21 (GeV 2 ) for COMPASS. One can see from Table UTl that description of the individual 
subsets is also very satisfactory. 

As we mentioned in Sec. [T] it is very important to include appropriate QCD evolution in order to be able to have 
a controlled accuracy and adequate description of e + e _ data that is measured at Q 2 ~ 110 GeV 2 . One can see from 
Table Hill th at both BELL E ITT and BABAR |98| data sets on Aq are described well. Both methods UC and UL from 
BELLE [12| and BABAR [98| appear to be consistent with our description and also Ph± dependence of asymmetry is 
well described. We will give predictions for Ph± dependent unpolarized cross sections in e+e - in the following section. 


Experiment 

hadron Target dependence # ndata 

x 2 

X 2 /ndata 

COMPASS [97] 

7T" 1 ” 

LiD 

X 

9 

11.16 

1.24 

COMPASS [97] 

7T _ 

LiD 

X 

9 

9.08 

1.01 

COMPASS [97] 

7T + 

LiD 

z 

8 

3.26 

0.41 

COMPASS [97] 

7 T~ 

LiD 

z 

8 

7.29 

0.91 

COMPASS [97] 

7T + 

LiD 

p h ± 

6 

4.19 

0.70 

COMPASS [97] 

7T _ 

LiD 

Ph± 

6 

4.50 

0.75 

COMPASS [96] 

7T -1 ” 

nh 3 

X 

9 

21.46 

2.38 

COMPASS [96] 

7V~ 

nh 3 

X 

9 

6.23 

0.69 

COMPASS [96] 

7V+ 

nh 3 

z 

8 

7.80 

0.98 

COMPASS [96] 

7 T~ 

nh 3 

z 

8 

10.29 

1.29 

COMPASS [96] 

7V + 

nh 3 

Ph± 

6 

3.82 

0.64 

COMPASS [96] 

7 T~ 

nh 3 

Ph± 

6 

3.85 

0.64 

HERMES [95] 

7V + 

H 

X 

7 

5.37 

0.77 

HERMES [95] 

7r - 

H 

X 

7 

12.61 

1.80 

HERMES [95] 

7T + 

H 

z 

7 

3.04 

0.43 

HERMES [95] 

7r _ 

H 

z 

7 

3.23 

0.46 

HERMES [95] 

7t + 

H 

p h ± 

6 

1.60 

0.27 

HERMES [95] 

7r _ 

H 

Ph± 

6 

4.82 

0.80 

JLAB [9] 

7r" 1 ” 

d He 

X 

4 

3.90 

0.98 

JLAB [9] 

7T _ 

3 He 

X 

4 

3.11 

0.78 





140 

130.65 

0.93 

TABLE II. Partial y 2 

values of the global best fit for SIDIS experiments. 


Experiment 

Observable dependence # ndata 

x 2 

y 2 /ndata 

BELLE [12] 

a ul 

Z 

16 

13.02 

0.81 

BELLE [12] 

aV c 

z 

16 

11.54 

0.72 

BABAR[ 98] 

aV l 

z 

36 

34.61 

0.96 

BABAR[ 98] 

A uc 

-A-0 

z 

36 

15.17 

0.42 

BABAR[ 98] 

A u L 

Ph± 

9 

9.09 

1.01 

BABAR[ 98] 

A uc 

Ph± 

9 

4.33 

0.48 


122 87.76 0.72 


TABLE III. Partial y 2 values of the global best fit for e + e experiments. 
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FIG. 3. Extracted transversity distribution (a) and Collins regimentation function (b) at three different scales Q 2 = 2.4 (dotted 
lines), Q 2 = 10 (solid lines) and Q 2 = 1000 (dashed lines) GeV 2 . The shaded region corresponds to our estimate of 90% C.L. 
error band at Q 2 = 10 GeV 2 . 



FIG. 4. x 2 profiles for up and down quark contributions to the tensor charge. The errors of points correspond to 90% C.L. 
interval at Q 2 = 10 GeV 2 . 


E. Transversity, Collins fragmentation functions and tensor charge 

We plot transversity and the Collins fragmentation function in Fig. [3] at two different scales Q 2 = 10 and 1000 
GeV 2 . In order to evaluate functions we solve appropriate DGLAP equations for transversity Eq. © and twist-3 
collins functions Eq. m- Due to the fact that neither of the functions mixes with gluons, these distributions do not 
change drastically in low-x region due to DGLAP evolution. 

Transversity enters directly in SIDIS asymmetry and we find that the main constraints come from SIDIS data only, 
its correlations with errors of Collins FF turn out to be numerically negligible. We thus vary only Xsidis an d use 
A Xsidis ~ 22.2 for 90% C.L. and A Xsidis = 6-4 for 68% C.L. calculated using Eq. (11231) . Since the experimental 
data has only probed the limited region 0.0065 < xb < 0.35, we define the following partial contribution to the tensor 
charge 


Sq [a 


min j^max 


/‘2-max 

1 (Q 2 ) = / dxh\(x , Q 2 ) . 

fCmin 


(127) 


In Fig. QJ we plot the x 2 Monte Carlo scanning of SIDIS data for the contribution to the tensor charge from such a 
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FIG. 5. x 2 profiles for up and down quark contributions to the tensor charge. The errors of points correspond to 68% C.L. 
interval at Q 2 = 10 GeV 2 . 


region, and find 


^[0.0065,0.36] = +0.30+0;“ } (128 ) 

^[0.0065,0.35] = -0.20™ , (129) 

at 90% C.L. at Q 2 = 10 GeV 2 . Analogously in Fig. [5j we plot the v 2 Monte Carlo scanning of SIDIS data at 68% 
C.L. at Q 2 = 10 GeV 2 and find 


^[0.0065,0.35] = +Q.3Q+0-04 j ( 130 ) 

^[0.0065,0.35] = _o.20t°;J 2 . (131) 


We notice t hat this result is comparable with previous TMD extractions without evolution PS® and di-hadron 
method [dl. Il03 ij. 

Existing experimental data covers a limited kinematic region, thus a simple extension of our fitted parametrization 
to the whole range of 0 < xb < 1 will significantly underestimate the uncertainties, in particular, in the dominant 
large-era regime. It is extremely important to extend the experimental study of the quark transversity distribution 
to both large and small Xb to constrain the total tensor ch arge contributions. This requires future expe riments t o 
provide measurements at the Jefferson Lab 12 GeV upgrade [l07j and the planned Electron Ion Collider I4S. ilOSL Il09l |. 
Nevertheless for completeness let us present our results on the tensor charge calculated over the whole kinematical 
region Jgl 0,1 !: 


W 0 ’ 1] = +0.39t°i® , 
Sd [ °’ 1] = -0.22±°;?J , 

at 90% C.L. and 

W°’ 1] = +0.39lg'll , 

<W [M = -0.22t“;St , 


(132) 

(133) 


(134) 

(135) 


at 68% C.L. both at Q 2 = 10 GeV 2 , as shown in Figs. [6] and [7) The tensor charge for u-quark can have a bigger 
contribution with respect to d-quark from unexplored region of x according to our estimates. 


F. TMD Interpretations of Our Results 


As we mentioned in the Introduction and elaborated in Sec. II, there exists TMD interpretation of CSS results. In 
particular the equations of the previous section that represent the solution of evolution equations are quite complicated. 
One might formulate the solutions of TMD evolution equations for TMD functions directly, in su ch a way that 
expressions will look very much like extension of a simple parton-like model, for instance used in Ref. m. 
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FIG. 6. x 2 profiles for up and down quark contributions to 
points correspond to 90% C.L. interval. 



tensor charge in the whole kinematical region. The errors of 




FIG. 7. x 2 profiles for up and down quark contributions to the tensor charge in the whole kinematical region. The errors of 
points correspond to 68% C.L. interval. 


Let us start with writing the cross sections in terms of the individual TMDs: 


Fuu{Q\ Ph±) 

FuT Ma \Q-,P h ±) 

Zuu h2 {Q\Ph±) 

^coUins (Q;Ph±) 


= HsimsiQ, H = Q) 1 


'k± ,p _l 


Q)D h/q (z h ,p\;Q), 


= —HsimsiQ, P> = Q) e a [ K{ x B,k±-,Q)-^j-r^H^ h / q (zh,p 2 j_-,Q), 

„ Jk±,p± M-h 


— He+e~ (Qi ft — Q) ^ ^ ( 


Dh/Q)Dh/q(Zh2i P2A-i 


(Q, n = Q) E e l [ ( 2p h± p L - 9?) 

„ J Pi I ,X?9 I X 7 


Pl± H lh/o( z hl,P 2 l±-,Q) 


'Pl±,P2± 
/2 


Mhi 


X J^ H lh/q( Z h2,p 2 2 J_]Q), 


(136) 


where we have set the factorization scale p = Q, and the evolution effects have been fully taken into account in the 
TMDs, fi(xB,kj_',Q), h\{x B ,k\\Q ), D h / q {z,p\\Q), and H^~ h ,-(z,p\;Q) are the transverse momentum dependent 
unpolarized parton distribution function, quark transversity, unpolarized fragmentation function, and the Collins 
function at the scale p = Q and C, = Q 2 , respectively. These TMDs are also understood as the soft factor subtracted 
TMDs. 

The short-hand notations for the integrations have the following explicit forms: 


J = J d 2 k±d 2 p±8 2 [zhk_L +p±- 


= / d 2 pi±d 2 p 2 ±8 2 ( P h _L ~Pl± ~P2± — 

Zh 2 


(137) 

(138) 


As discussed in Sec. II, the TMDs and the associated hard factors depend on the scheme to regulate the light-cone 
singularities. However, in the final results for the structure functions, this scheme dependence cancels out between 
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the TMDs and the hard factors. In the following, we present the results in the Collins-11 scheme [22}. The functions 
that encode scheme dependence from Eqs. (133I34I60I61I) are F q {a s (Q)) = 1, r D q (a s (Q)) = 1, Ji\ q {a s {Q)) = 1, 
' H c {a s {Q )) = 1 at one loop. For all other schemes, the results can be obtained accordingly. 

In the Collins-11 scheme, the associated hard factors can be written using Eqs. (1A20I98I) as 


#srois(Q, P = Q) = Id-(-8) , 

H e+e - (Q,n = Q) = 1 + C F (tt 2 - 8) . (139) 

The TMDs are Fourier transformations of the relevant expressions in b space in Sec. II, Eqs. (133I34I60I61I1 . 


, k 2 ± ;Q) = ^Mk±b) C f q U ® f[{x, p b ) e - 5 

h\{x,k\]Q) = f ^-J 0 (k±b) SC q ^ ^ h\(x, ii b ) e~^ Spert{Q ’ bt ' , ~ s ^p {Q ’ b \ 

Jo [zn) 

D h/q {z,pl-Q)= ^ [ ^Jo(p±b/z) C°1 ® D h/i (z, p b ) , 

z Jo \ Z7T ) 

j^ H th /q ( z ’Pl; Q) = 72 J o J^Mp±W z ) SC^ ns « fl$(z, p b ) e-^pe rt (Q,6,)-s^ lms (Q,6) ; 


(140) 

(141) 

(142) 

(143) 


where the TMD evolution has been taken into account, and one-loop results of T q , V q , "Hi q , and T~L C equal to one 
in the Collins-11 scheme have been applied, and C-functions are given in Eqs. (139I40I41I42I62I631) . Using relation to 
Trento conventions of Eq. (15811 we can write 


g-ffi/.O.pi; m—jsf t6g? « , 

(144) 

db b Ji{p±b/z ) <g> H^](z)\Trento( z , Pb) . 

(145) 


Hlh/q( Z iP-L’Q)\T' rent o — 


(2tt) 


We can also write explicitly the non-perturbative Sudakov form factor <Snp(Q, b) for all the TMDs discussed in our 
paper: 


S&(Q, b ) = Sxp(Q, b) = | In In + g q b 2 , 

In (£)+£*, 


(146) 

(147) 

(148) 


where we have assumed that the non-perturbative Sudakov form factors are the same for ff and h\ as a first study 
following With the expressions for individual TMDs given in Eqs. (114011141111421 fl45l) . and the fitted parameters 
in this section, we are now ready to present all these TMDs as a function of both longitudinal momentum fraction (x 
or z) and the transverse component ( k± or p±). 

In Fig. [5] we present unpolarized w-quark distribution /i at x = 0.1 as a function of b (left) and k± (right). We plot 


fi(x,b;Q ) 


b r h 

(2t r) ‘ 


fi(x,gb) e - 5 s 'pert(Q,b.)-5 , Np(Q- b ) ) 


(149) 


while fcj_-dependence is defined in Eq. (I140D . The distribution is calculated at three different scales: Q 2 = 2.4 (dotted 
lines), Q 2 = 10 (solid lines) and Q 2 = 1000 (dashed lines) GeV 2 . As one can see, at large scale Q 2 = 1000 GeV 2 , the 
distribution is highly dominated by perturbative region of b < b max while at lower scales Q 2 = 2.4 and 10 GeV 2 the 
distribution is shifted towards large values of b ~ 2 4- 3 GeV -1 , in this region of b one needs to carefully account for 
non-perturbative effects of TMD evolution and intrinsic motion of quarks. The distribution in k± space is becoming 
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(b) 


FIG. 8. Unpolarised w-quark distribution as function of b (a) and as function of k± (b) at three different scales Q 2 = 2.4 
(dotted lines), Q 2 = 10 (solid lines) and Q 2 = 1000 (dashed lines) GeV 2 . 


wider with growth of Q 2 and has developed perturbative tail, while at low valu es of Q 2 it resembles gaussian-type 
parametrization used in tree level extractions, for instance that of Refs. mini. 

The same observation is true for transversity distribution. We present transversity u-quark distribution hi at 
x = 0.1 as function of b and k± in Fig. [9] We plot 

hi(x,b;Q) = 7 ^6C q ^ i ®h\(x, t ,b)e-i s ^ Q ’ b ^- s ^ Q ’ b \ (150) 

(2tt) 


while k± distribution is defined in Eq. (11411) . Note that coefficient functions for transversity distribution SC qi -i are 
different from those of unpolarized distribution. This difference affects the shape of distributions in b and kj_ space. 
Moreover the width of transversity can be different from that of unpolarized distribution as well, however features of 
TMD evolution are very similar in both cases. Generic results on transversity TMD evolution were also presented in 

Ref. H. 

Unpolarized fragmentation TMD as function of b is defined as 


Dh/q(Z-> Q) — 


— — —C Dl < 
z 2 ( 2n) l ^ q 


'D h/i (z,p b )e 


-§S'pert(Q,f>,)-'S£p(Q>U 


(151) 


and as function of p± it can be calculated using Eq. (11421) . In Fig. QlJlwe present unpolarized TMD FF at z = 0.4 
and at three different scales Q 2 = 2.4 (dotted lines), Q 2 = 10 (solid lines) and Q 2 = 1000 (dashed lines) GeV 2 . Again 
as in case of other TMDs above one observes widening of distributions in p± and shift towards lower values b of the 
maximum of the distribution with increase of Q 2 . In relatively low Q 2 region the effects of TMD evolution are quite 
moderate. 

Collins fragmentation function with evolution is presented for the first time in this paper. The b dependent function 
can be defined as 


H 1 h/q( Z i Q) I Trento — 


collins , 


1 JLsc 

z 2 (2tt) 


ITrento^, Hb) e ! S P-t(Q,G)-S^ lms (Q,b) 


(152) 


and p± dependent function is in Eq. (11451) . In Fig. [TT]we present TMD Collins FF at z = 0.4 and at three different 
scales Q 2 = 2.4 (dotted lines), Q 2 = 10 (solid lines) and Q 2 = 1000 (dashed lines) GeV 2 . One observes widening of 
distributions in p± and shift towards lower values b of the maximum of the distribution with increase of Q 2 . Note 
that TMD Collins FF has a kinematical zero due to the pre-factor p±/zMh- 


G. Description of the experimental data 

The description of the HERMES data [9^] is shown in Fig. [12] One can see that the description is good for all 
Xb, z, and Ph± dependencies. The formalism that we use is appropriate in the region of low Ph± and we limit our 
description by Ph _l < 0.8 GeV. The data is in the region of 1 < ( Q 2 ) < 6 (GeV 2 ). The estimate of the error band is 
presented as the shaded region. 
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FIG. 9. Transversity M-quark distribution as function of b (a) and as function of k± (b) at three different scales Q 2 = 2.4 
(dotted lines), Q 2 = 10 (solid lines) and Q 2 = 1000 (dashed lines) GeV 2 . 



(b) 


FIG. 10. Unpolarised FF u — > -k + as function of b (a) and as function of p± (b) at three different scales Q 2 = 2.4 (dotted lines), 
Q 2 = 10 (solid lines) and Q 2 = 1000 (dashed lines) GeV 2 . 


One can see from Fig. [12] that both data and the model obey kinematical suppression of asymmetries at low 
Zh, and Ph_ l- Additionally the data indicates that asymmetry becomes smaller in the region of small-a;^ and thus 
transversity becomes small in the small -xg region as well as can be seen in Fig. [3] (a). Positive asymmetry of n + 
production implies that the product of u-quark transversity and the favored Collins fragmentation function is positive. 
We choose the solution with positive M-quark transversity (the same sign as M-quark helicity distribution) and obtain 
favored Collins fragmentation function is positive, see Fig. [3] (b). Large negative asymmetry of ir~ production indicates 
that the so-called unfavored Collins fragmentation function is large and negative and indeed it is the case, see Fig. [3] 
(b). Measurements on proton targets are dominated by w-quark functions as far as e 2 /e 2 = 4, thus we have better 
precision for the extraction of u-quark transversity and tensor charge Su. 

The COMPASS data [96j,:97j| extend the region of resolution scale by a factor of three, ( Q 2 ) < 21 (GeV 2 ). We present 
results of our description in Fig. [13] Again we exclude the region of Ph± > 0.8 GeV where relation Ph±/{z) < Q is 
not satisfied. The COMPASS data extends the region of xb up to xb ~ 10 -2 and the measured asymmetry indicates 
that transversity is rather small in the small -x region. Indeed the extracted transversity shown in Fig. [3] (a) becomes 
small in the small-a; region. The COMPASS data on effective deuterium target Fig. [U2 (b) indicate that the sum of 
w-quark and d-quark transversities is small, and thus both functions are approximately of the same size, it can be 
seen in Fig. [3] (a). 

Description of JLab’s HALL A data (9j is shown in Fig. [Li] The data extend the region of Xb toward large-x 
and one can see that our fit is compatible with the data. The measurement on effective neutron target ( 3 He) is 
sensitive to d-quark functions, however the current experimental errors are too big to allow better extraction of 
d-quark transversity. 

Both BELLE [l2j and BABAR [9S| collaborations measured the Collins asymmetries in e + e _ at yfs — 10.6 GeV. 
Comparison of BELLE data [l2j on Aq asymmetries for both UL and UC methods is presented in Fig. [15] The data 
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FIG. 11. Collins FF u —> n + as function of b (a) and as function of p± (b) at three different scales Q 2 = 2.4 (dotted lines), 
Q 2 = 10 (solid lines) and Q 2 = 1000 (dashed lines) GeV 2 . 



FIG. 12. Description of Collins asymmetries measured by the HERMES Collaboration [95| as a function of xb, z, Ph± on the 
proton target. The shaded region corresponds to our estimate of 90% C.L. error band. 


are measured in four different bins of Zhi,Zh 2 and one can see that the description of the data is very good. The 
asymmetry becomes small when Zhi,Zh 2 become small due to kinematical suppression and one can see from Fig. 1151 
that our calculations are compatible with this behavior. 

In Fig. [16] we present description of BABAR data [98[ on Aq asymmetries for both UL and UC methods. The 
data are in six bins of Zhi , Zh 2 with six points in each bin. This allows for better extraction of the shape of Collins 
fragmentation functions. One can see that also in this case the description is very good. The large -2 region deserve 
a special comment. One expects that the formalism will become unreliable when Zhi —> 1 and/or Zh 2 —> 1 due to the 
influence of exclusive pion production. Indeed one can see from Figs. [Tbl and ITTjl that in large -2 bins the quality of 
description deteriorates. Nevertheless both magnitude and the shape of the data are reproduced perfectly in the plot. 
It is achieved by allowing parameters that describe shape of favored and unfavored Collins fragmentation functions 
be different and independent of each other. Additionally the correct Q 2 evolution reproduces the shape much better 
compared to the case of absence of the evolution. Note that we have not attempted to fit the data without TMD 
evolution, thus our conclusion is valid only for comparison of results with and without evolution using parameters of 
NLL fit. 

Even though a priori it is very difficult to expect perfect description of the data in the whole 2 region, our fit indeed 
is capable of reproducing the data very well. Both Aq L and Aq C are described very well, we observe no tension 
between the measurements and it indicates the robustness of the method. Aq L and Aq C have slightly different 
sensitivity to different combinations of Collins fragmentation functions as can be seen from Eq. ([8]) and the usage of 
both measurements helps to constrain the functions better. We believe that favored Collins fragmentation functions 
are well determined and future experimental data could test our findings. 

Finally we present comparison of our calculations with Phi. dependence of e + e _ asymmetries in Fig. 1171 Both A^ L 
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FIG. 13. Description of Collins asymmetries measured by the COMPASS Collaboration as a function of xb, z, Ph± on (a) NH3 
proton f9l| and (b) LiD deuterium |97j targets. The shaded region corresponds to our estimate of 90% C.L. error band. 



FIG. 14. Description of Collins asymmetries measured by the JLab’s HALL A [ 9 ] as a function of xb on the effective neutron 
target. The shaded region corresponds to our estimate of 90% C.L. error band. 


and A% G are described very well. 

From the comparisons of the data and theoretical computations we can deduce that the TMD evolution at NLL' 
can describe both e + e _ and SIDIS data adequately well. The highest resolution scale in our analysis is quite big 
Q 1 2 = s ~ 110 GeV 2 and we found that using appropriate QCD evolution was essential in order to describe the data. 
It allows us to have a controlled theoretical precision of our computations. Let us study sensitivity of our results to 
the theoretical precision of computations. We will fix the parameters to the NLL' fit results presented in Table. U 
and calculate asymmetries in different kinematical configurations using tree level approximation, i.e. without TMD 
evolution, Leading Logarithmic accuracy (LL), and Next-to-Leading Logarithmic accuracy NLL'. As far as parameters 
are defined by fitting at NLL' we expect that NLL' will describe the data better than LL or tree approximation. We 
will not attempt to fit data at either tree approximation or LL, even though such fits can be well performed and may 
give reasonable descriptions of the data. By computing results with three different precisions with fixed parameters 
we will be able to answer two different questions: 

1. How big are effects of inclusion of higher orders in calculation of a particular asymmetry in a particular kine¬ 
matical region? 

2. How sensitive are experimental data to the inclusion of higher orders? 

We show NLL', LL, and no TMD evolution results for asymmetry as function of xb, z, Ph± for HERMES in 
Fig. [TU The computation at Leading Logarithmic accuracy (LL) is done by using only A W in the perturbative 
Sudakov form factor and C ^ coefficient function. No TMD evolution implies that the perturbative Sudakov form 
factor and parameter (72 are set to zero, accordingly we use DSS LO for fragmentation functions and CTEQ6LO for 
distribution functions and set the scale to Q q = 2.4 GeV 2 . The dotted line in Fig. IT51 shows result without TMD 
evolution. One can see that at low energy, the results are quite similar for all three calculations. This happens 
due to the fact that in ratios most of the numerical effects of evolution cancel out. The precision of existing SIDIS 
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FIG. 15. Collins asymmetries in e + e _ at y/s = 10.6 GeV measured by BELLE collaboration [l3| as a function of Zh 2 in different 
bins of Zhi (a) UL and (b) UC. Calculations are performed with parameters from Table Q] The shaded region corresponds to 
our estimate of 90% C.L. error band. 
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FIG. 16. Collins asymmetries in e + e~ at yfs = 10.6 GeV measured by BABAR Collaboration j9j|] as a function of Zh 2 in different 
bins of zhi (a) UL and (b) UC. Calculations are performed with parameters from Table H) The shaded region corresponds to 
our estimate of 90% C.L. error band. 


experimental data is such that it does not allow to distinguish among different theoretical accuracies used to calculate 
TMDs. It happens due to the fact that both energy and Q 2 are quite low for SIDIS. The difference grows as we 
consider COMPASS data in Fig [HI and the data becomes sensitive to the choice of accuracy. One can also see from 
Fig. [IS] that the difference in different precisions (no evolution, LL, NLL’) is comparable with the error band of the 
NLL’ extraction. One can conclude that results of phenomenological extraction using different precisions will be very 
similar if low energy experiments are used. In fact low energy experimental data are dominated by non-perturbative 
physics and even tree level approximation (no TMD evolution) grasps well the features of the underlying physics. We 
will in fact see that our results compare very well to the results of Torino-Cagliari-JLab group Ha¬ 
lil Fig. [20] we show theoretical computations for e + e _ without TMD evolution (dotted line), LL accuracy (dashed 
line), and the complete NLL' accuracy (solid line). The difference between these computations diminishes when we 
include higher orders, it means that the theoretical uncertainty improves. We conjecture that the difference between 
NLL' and NNLL will be smaller than difference between NLL' and LL and thus be comparable to experimental errors. 
One can also observe that asymmetry at Q 2 = 110 GeV 2 is suppressed by factor 2-3 with respect to tree-level 
calculations due to the Sudakov form factor. One can also conclude that NLL accuracy is essential for e + e _ data. 
Notice that we present calculations with fixed parameters determined by the NLL' fit. The difference between different 
curves shows sensitivity to the theoretical accuracy and to inclusion of higher order. The observation that calculation 
without TMD evolution or LL cannot describe the data with these parameters does not mean that a fit of the data 
without TMD evolution or LL is impossible. In fact such fits are most probably possible and could yield results of 
similar quality of description of the data. There is no doubt however that higher order computations such as NLL have 
an advantage of having better control of theoretical uncertainty. The fact that we utilize NLO collinear distributions 
is very encouraging, these distributions describe inclusive data sets much better than LO distributions. We also 
observe that e + e _ experiments are very sensitive to the inclusion of higher order corrections. This can be clearly seen 
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FIG. 17. Collins asymmetries in e + e _ at y/s = 10.6 GeV measured by BABAR collaboration [98j] as a function of Ph± (a) UL 
and (b) UC. Calculations are performed with parameters from Table U) The shaded region corresponds to our estimate of 90% 
C.L. error band. 



FIG. 18. Collins asymmetries measured by HERMES Collaboration [95| as a function of xb, z, Ph± in production of 7r + (left 
panels) and n~ (right panels). The solid line corresponds to the full NLL' calculation, the dashed line to the LL calculation, 
and the dotted to the calculation without TMD evolution. Calculations are performed with parameters from Table [Q The 
shaded region corresponds to our estimate of 90% C.L. error band. 


from Fig. l2ll where we compute Collins asymmetries measured by BABAR [98j Collaboration as a function of Zh 2 in 
different bins of z^i- One can see that importance of higher orders increases with increasing value of z^- 


H. Predictions for unpolarized multiplicities in SIDIS and e + e 


We predict unpolarized cross section of charged pion production to be measured by BELLE, BABAR and BESIII 
Collaborations and given by the formula 


dcr o 

dPL 


1 ( ^3 ( je + e ~^>h 1 h 2 +X 

(1 + cos 2 9) dzhidzh 2 dP%j_ 


N cTT 2 al m 

2 Q 2 


7hill2 


(153) 


where hi, /12 can be any charged pion, Z\ = Z 2 = 0.3. The prediction is given in Fig. 1221 As one can see we predict that 
measured cross-section will be wider for BELLE and BABAR Collaborations Q 2 = 110 GeV 2 with respect to BESIII 
Collaboration Q 2 = 13 GeV 2 . At the same time BESIII Collaboration Q 2 = 13 GeV 2 cross section will be larger 
than that measured by BELLE and BABAR Collaborations at Q 2 = 110 GeV 2 . In Fig. [25] we divide the predicted 
cross-section for BESIII Collaboration by a factor 110 in order to compare widths with expected cords-section at 
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FIG. 19. Description of Collins asymmetries measured by the COMPASS Collaboration on NH 3 proton [9f| as a function of 
xb, z, Ph± in production of 7 r + (left panels) and n~ (right panels). The solid line corresponds to the full NLL' calculation, the 
dashed line to the LL calculation, and the dotted to the calculation without TMD evolution. Calculations are performed with 
parameters from Table H) The shaded region corresponds to our estimate of 90% C.L. error band. 
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FIG. 20. Collins asymmetries measured by BABAR [& 8 j collaboration as a function of Ph± in production of unlike sign “U” 
over like sign “L” (a) and charged “C” (b) pion pairs at Q 2 — 110 GeV 2 . The solid line corresponds to the full NLL' calculation, 
the dashed line to the LL calculation, and the dotted to the calculation without TMD evolution. Calculations are performed 
with parameters from Table Q] 


BELLE and BABAR Collaborations. Effective widening of the cross-section with growth of Q 2 is a sign of TMD 
evolution and the future data from BELLE and BABAR Collaborations will be crucial for our understanding of the 
evolution. 

Similar behavior is shown in Fig. [23] of the unpolarized cross-section predicted for EIC at ^fs = 70 GeV and at 
Q 2 = 10 GeV 2 and Q 2 = 100 GeV 2 , choosing (zh) = 0.36 and (y) = 0.53. We plot 

j4 

-— g 2 = TT<T 0 (x B ,y,Q 2 )F uu . (154) 

dxBdydzhd z Fh±- 

The ultimate test of the TMD evolution will be in measurements of unpolarized cross-sections. We highly encour¬ 
age BELLE, BABAR and BESIII Collaborations to perform the analysis of the data on unpolarized cross-sections. 
Such measurements will allow us to test predictions of TMD evolution and will allow for a better understanding of 
unpolarized TMD fragmentation functions that can be measured directly only at e + e _ facilities. 

The universality of TMD evolution will be further tested in future measurements at Electron Ion Collider. EIC 
can easily span several decades in Q 2 and allow for a much better understanding of the nucleon 3D structure. The 
data of EIC combined with that of Jefferson Lab 12 will cover a very wide region of x and provide a multi-binning 
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FIG. 21. Collins asymmetries measured by BABAR Collaboration as a function of Zhi in production of unlike sign “U” 
over like sign “L” pion pairs at Q 2 = 110 GeV 2 . The solid line corresponds to the full NLL' calculation, the dashed line to the 
LL calculation, and the dotted to the calculation without TMD evolution. Calculations are performed with parameters from 
TableU The shaded region corresponds to our estimate of 90% C.L. error band. 



FIG. 22. Prediction for unpolarized cross-section in e + e at Q 2 = 110 GeV 2 to be measured by BELLE and BABAR 
Collaborations (solid line) and BESIII Collaboration at Q 2 = 13 GeV 2 (divided by a factor 110, dashed line) as a function of 

p2 

Ci_L- 


data needed for future phenomenological analysis. We plan to study the impact of EIC and Jefferson Lab 12 data in 
future publications. 


I. Predictions for future experiments in SIDIS and e + e 

BESIII is collecting data EH in e + e at Q 2 ~ 13 GeV 2 . The preliminary rcults are compatible with bigger 
asymmetries predicted by two of us in Ref. [40]. Here we present updated predictions assuming the same binning as 
BABAR and the following values of (sin 2 0}/(l + cos 2 9) = 0.65 at each bin, we also integrate the result in the following 
region of Phi. < 1 GeV. Actual values of asymmetry will depend on details of binning and kinematics. The predictions 
are presented in Fig. 1241 We give predictions for Aq L asymmetries, we predict enhancement of the asymmetry by 
factor 2 — 3, compare to Fig. [TG] (a). Note that our predictions will have to be scaled with actual experimental values 
of (sin 2 9)/(1 + cos 2 9) exp from BESIII. 

Measurements of Collins asy mmetries a re going to be performed at Jefferson Lab 12 GeV upgrade [l07| and the 
planned Electron Ion Collider fi flM Ilfp. The high precision of Jefferson Lab 12 measurements will eventually allow 
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FIG. 23. Prediction for unpolarized cross-section for an Electron-Ion Collider of energy y/s = 70 GeV at Q 2 = 100 GeV 2 (solid 
line) and Q 2 = 10 GeV 2 (divided by a factor 200, dashed line) as a function of P 2 ±. 
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FIG. 24. Predictions for UL Collins asymmetries in e + e at Q 2 = 13 GeV 2 to be measured by BESIII [l 1 ll ] as a function of 
Zh 2 in different bins of Zhi- 


for better determination of transversity distributions in high-x region and low-x region along with higher span in Q 2 
will be covered by EIC. 

Electron Ion Collider is going to allow studies of evolution in Q 2 and energy y/s od single spin asymmetries. It is 
going to provide a big leverage arm in Q 2 and will have variable center of mass energy y/s. We present here predictions 
of Collins asymmetry as function of Xb for two different values of Q 2 = 10 GeV 2 and Q 2 = 100 GeV 2 in Fig. [25] Note 
that Xb and Q 2 are correlated via Q 2 = sxbU , we also fix average values of Zh and P^_l, ( Zh) = 0.36 and ( Ph±) = 0.4 
GeV. One can see from Fig. [55] that we predict a moderate decrease of the asymmetry with Q 2 . Measurements in 
low-x region are going to provide information on sea quark transversity. Our current extraction neglects sea quarks, 
so the asymmetry becomes very small in low-x region. 

The Jefferson Lab 12 GeV program is going to extend our knowledge of the underlying distributions in the large-x 
region. Both proton and neutron ( 3 He) targets will provide information of distributions of u and d-quarks. We present 
predictions for JLab 12 at 11 GeV incident electron beam on proton and 3 He (effective neutron) targets in Fig. [5B] 
One can see that we predict sizable asymmetries of order of 10%, future data is going to highly improve the knowledge 
of transversity in large-x region, currently the error band is very big, see Fig. 1261 In order to give predictions in Fig. 1261 
we fixed the average kinematical variables, ( y) = 0.57, ( Zh ) = 0.5, and ( Phj_) = 0.38 GeV. 


J. Comparison to other extractions 

Tree level extraction of transversity and Collins fragmentation functions was performed by Torino-Cagliari-JLab 
group in papers [15 Nt 7| . In Fig. [57] (a) we present comparison of extracted transversity at NLL and result of Ref. E3- 
We also compare to extraction of transversity via dihadron fragmentation method [18| Fig. [57] (b). One can see that 
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X B 

FIG. 25. Predictions for Collins asymmetry as a function of xb for 7r + production on protons in SIDIS at Q 2 = 10 GeV 2 (solid 
line and vertical-line hashed region) and Q 2 = 100 GeV 2 (dashed line and shaded region) to be measured by EIC at energy of 
y/s = 70 GeV. 



X B 


FIG. 26. Predictions for Collins asymmetry as a function of xb for 7r + production at Jefferson Lab 12 GeV on proton target 
(solid line and vertical-line hashed region) and effective neutron target (dashed line and shaded region). 


all three extractions give consistent results in the explored region of xb • Within error bands of each extraction results 
are compatible with each other. One can see that the experimental data indeed show some tension, saturation of 
Softer bound, for d-quark in high-a; region as predicted in Ref. (9lj| . This saturation happens in the region not explored 
by the current experimental data, so future data from Jefferson Lab 12 will be very important to test the Soffer bound 
and to constrain the transversity and tensor charge. 

The functions themselves are slightly different as can be seen by comparing solid and dashes lines in Fig. \2H 
(a). In fact Ref. pa uses tree level TMD expression (no TMD evolution) for extraction and we use NLL TMD 
formalism. Results should be different even though in asymmetries, as we saw, at low energies results with NLL 
TMD are comparable with tree level. At higher energies and Q 2 situation changes and extracted functions must be 
different. At the same time one should remember TMD evolution does not act as a universal Q 2 suppression factor. A 
complicated Fourier transform should be performed that mixes Q 2 and b dependence and thus the resulting functions 
are different in shape, but comparable in magnitude. It is also very encouraging that tree level TMD extractions 
yielded results very similar to our NLL extraction. This makes the previous phenomenological results valid even 
though the appropriate TMD evolution was not taken into account. It also means that we need to have experimental 
data on unpolarized cross sections differential in Phj_- As we have seen the effects of evolution should be evident in 
the data and those measurements will help to establish validity of the modern formulation of TMD evolution. 

We compare extracted Collins fragmentation functions —zH^ 3 \z) in Fig. [25] at Q 2 = 2.4 GeV 2 with extraction 
of Torino-Cagliari-JLab 2013 E3- The resulting Collins FF have the same signs but shapes and sizes are slightly 
different. Indeed one could expect it as far as Q 2 of e + e _ is different and evolution effect must be more evident. At 
the same time those functions for both tree level and NLL TMD give the same (or similar) theoretical asymmetries 
that are well compared to the experimental data of SIDIS and e + e _ . The favored Collins fragmentation function is 
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FIG. 27. (a) C omp arison of extracted transversity (solid lines and shaded region) Q 2 = 2.4 GeV 2 with Torino-Cagliari-JLab 

2013 extraction 47] (dashed lines and shaded region). 

(b) Comparison of extracted transversity (solid lines and shaded region) at Q 2 = 2.4 GeV 2 with Pavia 2015 extraction fl8ll 
(shaded region). 



FIG. 28. Comparison of extracted Collins fragmentation functions (solid lines) at Q 2 = 2.4 GeV 2 with Torino-Cagliari-JLab 
2013 extraction [l7] (dashed lines and shaded region). 


much better determined by the existing data, as one can see from Fig. [21] that the functions at Q 2 = 2.4 GeV 2 are 
compatible within error bands. The unfavored fragmentation functions are different, however those functions are not 
very well determined by existing experimental data. 

We also compare the tensor change from our and other extractions in Fig. [29] The contribution to tensor charge 
of Ref. [l8j is found by extraction using the so-called dihadron fragmentation function that couples to collinear 
transversity distribution. The corresponding functions have DGLAP type evolution known at LO and were used in 
Ref. [48j. The results plotted in Fig. [?9l corresponds to our estimates of the contribution to u-quark and d-quark in 
the region of x [0.065,0.35] at Q 2 = 10 GeV 2 at 68% C.L. (label 1) and the contribution to u-quark and d-quark in 
the same region of x and the same Q 2 using the so-called flexible scenario, a s (M|) = 0.125, of Ref. [18]. One can 
see that our extraction has an excellent precision for both it-quark and d-quark. The fact that the central values and 
errors of extracted tensor charges are in a good agreement in both methods, ours and Ref. [18], is very positive and 
allows for future investigations of transversity including all available data in a global fit. 

Our results compare well with extractions from Ref. G3- Even though correct TMD evolution was not used in 
Ref. 03 the effects of DGLAP evolution of collinear distributions were taken into account and the resulting fit is of 
good quality, x 2 /d.o./. = 0.8 for the so-called standard parametrization of Collins fragmentation functions. In fact 
the probability that the model of Ref. [TtJ correctly describes the data is P(0.8 * 249,249) = 99%. The tensor charge 
was estimated at 95% C.L. using two different parametrizations for Collins fragmentation functions, the so-called 
standard parametrization that utilized similar to our parametrization and the polynomial parametrization. In Fig. 1301 
we compare our results with calculations from Ref. [13] at 95% C.L. at Q 2 = 0.8 GeV 2 and calculations at 68 % at 
Q 2 = 1 GeV 2 of Ref. [l8[ . Even though we compare tensor charge at different values of Q 2 its evolution is quite slow, 
so the good agreement of all three methods is a good sign. We conclude that tensor charge perhaps is very stable with 
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FIG. 29. Comparison of tensor charge <5q'^ 0 ' 0065,0 ' 35 ^ for n-quark and d-quark from this paper at 68% C.L. (Kang et al 2015) 
and result from Ref. [lf| (Radici et al 2015) at 68% C.L. Both results are at Q 2 = 10 GeV 2 . 
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FIG. 30. Comparison of tensor charge for n-quark and d-quark in the whole region of x from this paper at 90% C.L. 

(Kang et al 2015) at Q 2 = 10 GeV 2 and result from Ref. flSf (Radici et al 2015) at at 68% C.L. and Q 2 = 1 GeV 2 , and Ref. [Til l 
at 95% C.L. standard and polynomial fit (Anselmino et al 2013) at Q 2 = 0.8 GeV 2 . 


respect to evolution effects that are included in phenomenological extractions. It also means that phenomenological 
results of Ref. EH and other extractions without TMD evolution are valid phenomenologically. One should remember, 
of course, that TMD evolution is more complicated if compared to DGLAP evolution (even though formal solutions 

are simpler in TMD case). The usage of non perturbative kernels make it very important to actually demonstrate 

that the proper evolution is indeed exhibited by the experimental data. Once correct evolution and non perturbative 
Sudakov factor are established the results of Ref. E2 should be improved by utilizing the appropriate TMD evolution 
that we have formulated in this paper. 

In Fig. [31] we compare tensor charge for u and d-q uarks from this paper at 90% C.L. at Q 2 = 10 GeV 2 

and resul ts from various model estimates of Refs. [nunf. One can see that our results are close to results of 
Ref. fl!3 j that actually used the approximate mass degeneracy of the light axial vector mesons (ai(1260), 1235) 

and /ii(1170)) and pole dominance to calculate the tensor charge. DSE calculations of tensor charge of Ref. [l!2l ] are 
also close to our results. 

Finally we present our estimates for the isovector nucleon tensor charge gr = Su — 5d: 

9t = +0.6lt°; 2 5 f , (155) 

at 90% C.L. and 

9 t = +0.6lto.25 . (156) 

at 68% C.L.at Q 2 = 10 GeV 2 . This result can be compared to lattice QCD calculations. 
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FIG. 31. Comparison of tensor charge for u-quark and rf-qu ark in the whole region of x from this paper at 90% C.L. 

(Kang et al 2015) at Q 2 — 10 GeV 2 and results from Refs. Ill2l4lld ]. 


In Fig. [321 we compare our result with extraction of Radici et al Ref. [l8} at Q 2 = 4 GeV 2 , Anselmi no e t al Ref. El 
standard and polynomial at Q 2 = 0.8 GeV 2 , and a series of lattice computations. Bali et al Re f. m estimate gx 
at TOtt ~ 150 MeV using RQCD with 2 flavor NPI Wilson-clover fermions, Gupta et al Ref. fm use 2 + 1 + 1 flavor 
HISQ lattices generated by the MILC collaboration with lowest = 130 MeV, Green et al Ref. El use 2 + 1 flavor 
BMW clover-improved Wilson action with pion masses between 149 and 356 MeV, Aoki et al use gauge configurations 
generated by the RBC and UKQCD Collaborations with (2 + l)-flavor QCD with domain wall fermions, PNDME 
Collaboration Bhattacharya at al |l20j use wo ensembles of highly improved staggered quarks lattices generated by 
the MILC collaboration with 2 + 1 + 1 dynamical flavors at a lattice spacing of 0.12 fm and with light- quark masses 
corresponding to pio ns w ith ma sses 310 and 220 MeV. references to other calculations of gx on lattice can be found 
for instance in Ref. |l20| . Ref. |l2l| uses rif = 2 lattice QCD, based on clover-improved Wilson fermions. One can 
see from Fig. 1321 that all phenomenological extractions indicate small values for the isovector nucleon tensor charge 
compared to lattice QCD. DSE computations of gx at Q 2 = 4 GeV 2 were performed in Ref. [l!2| and the result is 
different from most of lattice computations and closer to phenomenological extraction from the data. 

The value of gx extracted from the data may influence searches of BSM physics that depend on gx [l22l . Il23j . One 
can see that our determination of gx is the most precise existing extraction from experimental data. 

The isoscalar nucleon tensor charge g^ = Su + Sd can be readily computed using our results. We present result for 
gj, for completeness 

9t = +0.17±g;|S, (157) 

at 90% C.L. at_Q 2 = 10 GeV 2 . 

Refs, n+iin+f explores large- N c behavior of parton distributions in QCD and predicts that 

| hi(x) - hf(x )| » \hi(x) + hf(x )\, (158) 

we indeed observe that transversity for u and d-quarks are of similar magnitude and opposite signs and gx > gx and 
thus our results are compatible with large- N c predictions. 


IV. SUMMARY 

In this paper, we have performed a global analysis of the Collins azimuthal asymmetries in e + e _ annihilation and 
SIDIS processes, for the first time, with full QCD dynamics taken into account, including the appropriate TMD 
evolution effects at the NLL' order and perturbative QCD corrections at the NLO. The valence quark contributions 
to the nucleon tensor charge were estimated based on our analysis. Let us summarize the major results of this 
comprehensive study. 

First, the full QCD evolution effects are crucial to describe the Collins asymmetries in the back-to-back di-hadron 
productions in e + e _ annihilations, where current data come from the 13-factories at the center of mass energy around 
10.6 GeV. At this energy range, the TMD evolution has significant effect on the asymmetry distributions as functions 
of the transverse momentum, and the longitudinal momentum fractions carried by the hadrons in the fragmentation 
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FIG. 32. Comparison of the isovector nucleon tensor charge gr from this paper at 68% C.L. (Kang et al 2015) at Q 2 = 10 
GeV 2 and result from Ref. [l8j (Radici et al 2015) at 68% CL and Q 2 — 4 GeV 2 , and Ref. [13] at 95% CL standard and 
poly nomi al fit (Anselmino et al 2 013) at Q 2 = 0.8 GeV 2 . Other points are latt ice computation at Q 2 = 4 Ge V 2 of Bali et al 
Ref. IllTj, Gupta et al Ref. Ill8ll . Green et al Ref. [l 19t] . Aoki et al Ref. Il27tl . Bhattacharya et al ref. Il20ll . Gockeler et al 
Ref. 1 12 lfl . Pitschmann et al is DSE calculation at Q 2 = 4 GeV 2 Ref. [l!2ll . 


processes. These features have been clearly demonstrated in Figs. I20H211 In particular, the transverse momentum 
dependence illustrates the effects coming from the Sudakov resummation form factors where the perturbative part 
plays an important role due to large value of the resolution scale Q ~ 10.6 (GeV). The associated scale evolution 
effects in the H^ 3 \z) is another important aspect in the calculations. The evolution kernel is different from that of 
the unpolarized fragmentation function, and it changes the functional form dependence of Zhi and Zh 2 - In addition, 
there is cancellation between favored and unfavored Collins fragmentation functions, not only the shape but also the 
size are modified with the full evolution effects taken into account. 

Second, because of relative narrow Q 2 range in the current SIDIS data, the evolution effects are not so evident as 
compared to that in e + e _ annihilation processes. This was shown in Figs. [lH] and 1191 However, we would like to 
emphasize that, in order to precisely constrain the quark transversity distributions, we need to perform the complete 
QCD evolution in the theoretical calculations of the asymmetries to compare to the experimental data. This will 
become more important with high precisi on data fr om future experiments at the Jefferson Lab 12 GeV upgrade |l07l ] 
and the planned Electron Ion Collider (3, 11081 . Il09( |. 

Third, the quark transversity distributions from our analysis are comparable to previous determinations, including 
the leading order analysis of the same Collins asymmetries in SIDIS and e + e~ annihilation processes, and the di¬ 
hadron fragmentation channel in DIS and e + e _ processes, see Fig. [27] In particular, the consistency between the 
Collins asymmetry analysis and the di-hadron fragmentation analysis is a strong encouragement toward a future global 
fit to include all experimental data to constrain the quark transversity distributions. 

We observe, however, the Collins fragmentation functions from our analysis are quite different from those determined 
from the leading order analysis in Ref. although they are in the same order of magnitude. To further test the 
evolution effects, we emphasize the importance of future experiment measurements, in particular, in the energy range 
different from R-factories, such as those from the BEPC II at the experiment BESIII. We have made predictions for 
these experiments in Figs. [?51 and [Ml We hope the data will become available soon, and can be included into the 
global fit in the near future. We encourage BELLE, BABAR and BESIII Collaborations to perform the analysis of the 
data on unpolarized cross-sections as such data are curtail for our understanding of TMD fragmentation functions. 

Finally, we summarize the nucleon tensor charge contribution from our analysis, 

^[0.0065,0.35] = +0.30^;^ , (159) 

^[0.0065,0.35] = _o.20+0;28 , (160) 

at 90% C.L. at Q 2 = 10 GeV 2 , in the kinematic range covered by the current experiments. 

^[0 0065 . 0 - 35 ] = + 0 . 30+0 04 } 

^[0.0065,0.35] = _0.20±°$ , 


at 68% C.L. at Q 2 = 10 GeV 2 . 


(161) 

(162) 
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If we extend to the complete x range by assuming the model dependence in our fit, we would obtain 

<W°' 1] = +0.39±g;ao , (163) 

<W [tu l = -0.22±°;?J , (164) 

at 90% C.L. at Q 2 = 10 GeV 2 and 

<W°’ 1] = +0.39lo'ii , (165) 

8d [0 ' 1] = -0.22±g;J| , (166) 


at 68% C.L. at Q 2 = 10 GeV 2 . We emphasize that the above constraints depend on the functional form in our 
analysis, and the numbers quoted here should be taken cautiously. It is, nevertheless, interesting to compare to 
previous determinations. In Fig. [57] we show the comparisons of the nucleon tensor charges between our results and 
other determinations, together with some model calculations and the lattice computations. 

Much of improvements can be made in the future. First, more experimental data are in horizon from the 12 GeV 
upgrade of Jefferson Lab experiments, which actually will cover large-a; region and is of crucial importance to constrain 
the quark transversity distribution in that region. Since the nucleon tensor charge contribution is an integral of the 
quark transversity distribution, future Jefferson Lab data will be very important to reduce the uncertainties quoted 
above, and the uncertainties we can not address at the moment, such as the kinematic extension to obtain Jgl 0,1 !. 

The TMD evolution and the procedure to perform the global analysis will be an important part in the future 
analysis for other observables, for example, the Sivers asymmetries in SIDIS. We plan to carry out this analysis in a 
future publication. 

A number of improvements can be pursued in the theoretical part of the formalism. In this paper, we have taken 
the approximate evolution kernel for the twist-three quark-gluon-quark correlation contribution to the fragmentation 
function H^ 3 \z). For a complete analysis, we should include other terms in this evolution equation. Although it 
may not be possible to have a closed evolution equations for both H( 3 \z) and the related twist-three fragmentation 
functions Hd(zi, Z 2 ), one should be able to estimate the contributions from these additional terms. Second, with 
more experimental data available, we shall include the flavor dependence in the non-perturbative form factors in the 
Collins fragmentation function in the CSS resummation formalism. In this paper, we have assumed that they are 
flavor independent. The flavor dependence of distribution and fragmentation functions will be explored in the future 
analysis with more data available, in particular, the data on the transverse momentum dependence of the asymmetries 
in e + e^ annihilation processes. 

As a final remark, we would like to emphasize that our results and the methodology in the analysis will play 
an important role in phenomenological applications of perturbative QCD to the vast experimental data on SIDIS, 
Drell-Yan and e + e _ and in extraction of the relevant TMD parton distributions of the nucleon. 
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Appendix A: One-loop Calculation of the Collins Asymmetry in SIDIS 


To study the perturbative corrections and extract the hard factor in the Eq. (1551) , we need to carry out a calculation 
for E co ni ns at one-loop order. The leading order expression and the virtual diagram contributions follow that in the 
previous calculations for, e.g., the Sivers single spin asymmetry in SIDIS [39j, [52[. For the real gluon radiation, we use 
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the results in Ref. 
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where £ = Xb/x , £ = Zh/z and we only keep the most important diagonal contributions from H h j^(z). and the 

contributions from Hd(zi,Z 2 ) can be found from Ref. [HI|. By applying the Fourier transform (some of the useful 
integrals are listed in the Appendix of Ref. Q and Eq. (42) of Ref. [l28l p. we obtain the following result for 
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where we have partial integrated out the derivative terms in the previous equation to simplify the above expression. 
Clearly, the real diagrams contributions contain soft divergence (1/e 2 ), which will be cancelled by the virtual diagrams 
contributions. The virtual diagram contributes to a factor, 
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After canceling out these divergences, we have the total contribution at one-loop order, 
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where P represents the associated splitting kernels. They can be derived from the above results, 
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where we only list the part we have shown in the above from the contribution from H^ 3 \z) term. In general, the 
evolution of twist-three correlation functions involves multiple parton correlation contributions, for which there is no 
homogenous form. 
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To demonstrate the TMD factorization and calculate the hard factor in the TMD factorization, we have to calculate 
the transverse momentum dependence in the quark transversity distribution and the Collins fragmentation function 
at one-loop order. For the transversity distribution, we have 
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for the un-subtracted distribution in the JMY scheme. Adding the virtual contribution, 
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we obtain the total contribution for the un-subtracted quark transversity TMD at the one-loop order, 
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in the JMY scheme. Therefore, the subtracted TMD quark transversity distribution can be written as 
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(A7) 


(A8) 


(A9) 


,(A10) 


where we have also applied the renormalization for the integrated transversity distribution. By setting x B C, 2 = pp 2 
and p = p b as the initial scales for the TMD evolutions, we obtain the C-coefficient and the hard function Hi q of 
Eq. (l60l) as 


n[ J q MY) = i + ^c F 


27r 


lnp--ln p- — -2 


, 6C q ^ q (Z,p b ) = 6(l-Z)(l+0(a 2 s )) . (All) 


Similarly, we can carry out the calculations in the JCC scheme, for which we have the TMD quark transversity at 
one-loop order, 


7 SUb(JCC) ( 7 A \ 

/ii (x b , o; C? f 1 ) 




+c F d(i - 0 
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- - In 
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(A12) 


Applying the above result into Eq. (RjUI) and setting ( c = P = Pb as the initial scales for the TMD evolutions, we have 


H 


(JCC) 

D 


= l + 0(a 2 ), 6C q ^ q (Z,p b ) = 6(l-0(l + O(a 2 s )) . 


(A13) 


The above calculations can be extended to the Collins fragmentation function. The transverse momentum dependence 
can be calculated from perturbative QCD, and written in terms of the twist-three fragmentation function, 


Hl(Zh,P±) 


= p*c F . 1 


dz 


real 


2k 2 * (p 2 ± ) 2 J z 


-2f 6 


d H$ q {z) 
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(1-6+ v p 2 ± 
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Fourier transforming into 6 -space and adding the virtual diagram contribution (similar to that in Ea. (IA 8 l) l. we obtain 
the un-subtracted Collins fragmentation function at one-loop order, 
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2 c 2 p 2 2 


- 2 - 


in the JMY scheme. For the subtracted Collins fragmentation function, we have, 
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Similarly, we obtain the subtracted Collins fragmentation function in the JCC scheme, 
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From the above results, we derive the associated C-functions, 


(A15) 


• (A16) 


(A17) 


Ui JMY ) = 1 + ^ C F 

27r 


1 2 7T 2 

In p-In 2 p- 2 

H 2 2 


, p b ) = 5(1 - 0 + ^C F P^ q (0 Inf 


ni JCC > = 1 + 0 (a 2 ), 5 C ? ^(|, p 6 ) = 5(1 - i) + ^-C F P^ q (i) Inf . 

Finally, we can obtain the hard factors in both schemes. For example, in the Ji-Ma-Yuan scheme, 


^ mv (Q; m) = rtT )JMV (Q; m) = i + ^c f 


In + In p 2 In ^ — In p 2 + In 2 p + 7 r 2 — 4 

p 2 p 2 


(A18) 


(A19) 


Note that the hard part is the same for Fuu and Fjjt that is why we used the same notation FI in Eqs. (11715311 
Similarly, for the Collins-11 TMD scheme, we have 


*££i JOO (Q; aO = 4 D J S) JCC (Q; /*) = i + %c F 


3 In % - In 2 


Q 5 


(A20) 


These hard factors can be calculated from the factorization of E c “ llins , or from simply the virtual graphs for both the 
cross sections and the parton distribution and fragmentation functions. We will get the consistent results. 

In the end, the C-functions in Eq. (53) can be calculated from the above results, 


5C( SIDI S)(0 = 5C(f X n lq x V^coTnf^ = Q) 
5C( SIDI S)(|) = 5C(|) X H c X - Q) , 


(A21) 

(A22) 


where the scheme dependence is cancelled out between FLi q and In particular, the p dependence disappear 

in the JMY scheme when applying the above formulas to calculate the C-functions in the standard CSS resummation. 
Similarly, we can calculate the C-functions for the e"*"e _ annihilation processes, 


C (e+e_) (f = C(i) xV q x \JFduu e ~\p = Q) , 
5C'(e +e_) (f = 5C(f xil c x ^H^J(p = Q ) . 


(A23) 

(A24) 
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Again, the scheme dependence is cancelled out between the last two factors in the above equations. Comparing the 
SIDIS and e + e _ processes, we also find out that the difference comes from the hard factors. 
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